# ======================================================================================================== # # STEPHANE TUFFERY # DATA SCIENCE, STATISTIQUE ET MACHINE LEARNING - 6e édition # EDITIONS TECHNIP # 2025 # # CODE R ET PYTHON du livre # ======================================================================================================== # # --------------------------------------------------------------------------------------------------------- # Chapitre 1 : Data science et informatique # --------------------------------------------------------------------------------------------------------- # jeu de données Titanic install.packages('explore') library(explore) titanic <- as.data.frame(use_data_titanic(count = FALSE)) # remplacement des chaînes de caractère par des numéros de classe titanic[] <- lapply(titanic, function(x) as.integer(as.factor(x))-1) # on veut que la classe "crew" ("équipage") ait le numéro 0, la 1ère classe le numéro 1, etc. titanic$Class <- (1+titanic$Class)%%4 titanic$Age <- 1-titanic$Age head(titanic) # régression logistique titanic.model <- glm(Survived ~ . , data=titanic) coefficients(titanic.model) # export en PMML install.packages('pmml') library(pmml) pmml(titanic.model) # export en SQL library(tidypredict) tidypredict_sql(titanic.model, dbplyr::simulate_mssql()) # arbre de décision library(rpart) titanic.model <- rpart(Survived ~ . , data = titanic, method="class") print(titanic.model) # affichage de l'arbre avec package rattle library(rattle) fancyRpartPlot(titanic.model, sub=" ", type=1) library(rpart.plot) rpart.rules(titanic.model) # --------------------------------------------------------------------------------------------------------- # Chapitre 3 : Gestion des valeurs manquantes # --------------------------------------------------------------------------------------------------------- # package VIM install.packages('VIM') library(VIM) data(sleep, package = "VIM") # jeu de données "Sleep in Mammals" : http://www.statsci.org/data/general/sleep.html summary(sleep) # observations complètes sleep[complete.cases(sleep),] # nombre d’observations complètes sum(complete.cases(sleep)) # observations avec valeurs manquantes sleep[!complete.cases(sleep),] # visualisation des valeurs manquantes a <- aggr(sleep) a summary(a) aggr(sleep, prop = FALSE, numbers = TRUE) a <- aggr(sleep, prop = F, numbers = T, col=c("grey80", "grey60")) # graphique avec valeurs manquantes en rouge par défaut, et valeurs plus grandes en couleur plus foncéee matrixplot(sleep, interactive = FALSE, sortby = "Sleep") matrixplot(sleep, interactive = FALSE, sortby = "NonD", col="white") # graphique avec points rouges pour valeurs manquantes (aucune valeur manquante pour les deux variables simultanément) marginplot(sleep[,c("Gest","Dream")]) # imputation sleep_IMPUTED <- kNN(sleep) sleep_IMPUTED <- hotdeck(sleep) a <- aggr(sleep_IMPUTED, delimiter="_imp") a summary(a) vars <- c("Span","Dream","Dream_imp") marginplot(sleep_IMPUTED[,vars], delimiter="_imp", pch=c(19), col=c("grey70", "grey30")) marginplot(sleep_IMPUTED[,vars], delimiter="_imp", pch=c(19)) # --------------------------------------------------------------------------------------------------------- # Chapitre 3 : Distribution normale # --------------------------------------------------------------------------------------------------------- # Avec PYTHON : Figure 3.6 #----------------------------------------- import matplotlib.pyplot as plt import numpy as np # Generate data for the Gaussian curve x = np.linspace(-4, 4, 100) y = np.exp(-x**2 / 2) / np.sqrt(2 * np.pi) # Create the figure plt.figure(figsize=(10, 6)) # Plot the Gaussian curve plt.plot(x, y, color='black', linewidth=2) # Define the standard deviations and probabilities sigmas = [1, 2, 3] probabilities = [0.682, 0.954, 0.997] colors = ['#111111', 'gray', 'lightgray'] # Add vertical lines and shaded areas for sigma, probability, color in zip(sigmas, probabilities, colors): x_sigma = sigma plt.axvline(x=x_sigma, color='gray', linestyle=':') plt.axvline(x=-x_sigma, color='gray', linestyle=':') margin = 0.05 if sigma == 1 else 0 plt.fill_between(x, y, where=(x > -(x_sigma + margin)) & (x < (x_sigma + margin)), color=color, alpha=0.5) # Add text labels on the x-axis plt.xticks([-3, -2, -1, 1, 2, 3]) plt.text(-3, -0.037, "-3\u03C3", ha='center') plt.text(-2, -0.037, "-2\u03C3", ha='center') plt.text(-1, -0.037, "-1\u03C3", ha='center') plt.text(1, -0.037, "1\u03C3", ha='center') plt.text(2, -0.037, "2\u03C3", ha='center') plt.text(3, -0.037, "3\u03C3", ha='center') # Hide numerical x-axis ticks plt.xticks([-3, -2, -1, 1, 2, 3], ['' for _ in range(6)]) # Add horizontal line with horizontal arrows for 1 standard deviation plt.plot([-1, 1], [0.28, 0.28], color='black', linewidth=1) plt.quiver(-0.9, 0.28, -0.1, 0, angles = 'xy', scale_units = 'xy', scale = 1) # Tracé d'un vecteur plt.quiver(0.9, 0.28, 0.1, 0, angles = 'xy', scale_units = 'xy', scale = 1) # Tracé d'un vecteur # For 2 standard deviations plt.plot([-2, 2], [0.33, 0.33], color='black', linewidth=1) plt.quiver(-1.9, 0.33, -0.1, 0, angles = 'xy', scale_units = 'xy', scale = 1) plt.quiver(1.9, 0.33, 0.1, 0, angles = 'xy', scale_units = 'xy', scale = 1) # For 3 standard deviations plt.plot([-3, 3], [0.38, 0.38], color='black', linewidth=1) plt.quiver(-2.9, 0.38, -0.1, 0, angles = 'xy', scale_units = 'xy', scale = 1) plt.quiver(2.9, 0.38, 0.1, 0, angles = 'xy', scale_units = 'xy', scale = 1) # Add probability text plt.text(0, 0.4, f"99.7%", ha='center') plt.text(0, 0.35, f"95.4%", ha='center') plt.text(0, 0.3, f"68.2%", ha='center') # Adjust limits and labels plt.xlim(-4, 4) #plt.ylim(-0.1, 0.45) plt.ylim(-0.01, 0.45) plt.xlabel('0') plt.ylabel('Densité de probabilité') plt.yticks(np.arange(0, 0.45, 0.1)) # Save the plot in black and white plt.savefig('gauss_curve.png', bbox_inches='tight', dpi=300) # Télécharger la figure from google.colab import files files.download("gauss_curve.png") # Show the plot plt.show() # Avec R #----------------------------------------- plot(seq(-3.2,3.2,length=50), dnorm(seq(-3,3,length=50),0,1), type="l", xlab="", ylab="", ylim=c(0,0.5)) segments(x0 = c(-3,3), y0 = c(-1,-1), x1 = c(-3,3), y1=c(1,1), lty=2) text(x=0,y=0.45, labels = expression("99.7% within 3" ~ sigma)) arrows(x0=c(-1,1), y0=c(0.45,0.45), x1=c(-3,3), y1=c(0.45,0.45)) segments(x0 = c(-2,2), y0 = c(-1,-1), x1 = c(-2,2), y1=c(0.4,0.4), lty=2) text(x=0,y=0.3, labels = expression("95.4% within 2" ~ sigma)) arrows(x0=c(-1,1), y0=c(0.3,0.3), x1=c(-2,2), y1=c(0.3,0.3)) segments(x0 = c(-1,1), y0 = c(-1,-1), x1 = c(-1,1), y1=c(0.25,0.25), lty=2) text(x=0,y=0.15, labels = expression("68.2% within 1" * sigma)) arrows(x0=c(-0.65,0.65), y0=c(0.15,0.15), x1=c(-1,1), y1=c(0.15,0.15)) # --------------------------------------------------------------------------------------------------------- # Chapitre 7 : Analyse factorielle # --------------------------------------------------------------------------------------------------------- # lecture fichier SAS library(sas7bdat) etude <- read.sas7bdat("C:/Users/.../enquete.sas7bdat") library(haven) etude <- read_sas("C:/Users/.../enquete.sas7bdat") # lecture fichier CSV etude <- read.csv2("C:/Users/.../enquete.csv", stringsAsFactors = TRUE) head(etude) summary(etude) afc <- etude[,c("age","rayon")] summary(afc) table(afc$age, afc$rayon) # nombre de modalités par variable (modal <- apply(afc, 2, function(x) nlevels(as.factor(x)))) # nb de valeurs propres non triviales # = nb modalités - nb variables sum(modal) - ncol(afc) # inertie totale = (# modalités / # variables) - 1 (sum(modal)/ncol(afc))-1 #install.packages('RcmdrPlugin.FactoMineR') #library(RcmdrPlugin.FactoMineR) #library(Rcmdr) install.packages('FactoMineR') library(FactoMineR) # ACM # les variables doivent être catégorielles (ce n'est pas le cas avec l'import de SAS) afc$age <- as.factor(afc$age) afc$rayon <- as.factor(afc$rayon) # tableau de contingence multiple ACM <- MCA(afc, ncp=sum(modal), axes = c(1,2), graph = TRUE, method="Burt") # sur tableau disjonctif complet ACM <- MCA(afc, ncp=sum(modal), axes = c(1,2), graph = TRUE) AFC <- CA(table(afc$age, afc$rayon), graph = TRUE) AFC <- CA(table(afc$age, afc$rayon)) summary(AFC) # variables et modalités caractérisant chaque axe dimdesc(AFC) # test ANOVA pour mesurer la liaison entre chaque variable qualitative et chaque axe # valeurs propres AFC$eig barplot(AFC$eig[,2], names=paste("Dim",1:nrow(AFC$eig))) ACM$eig # inertie totale = somme des valeurs propres sur ACM sur le tableau disjonctif complet sum(ACM$eig[,1]) # somme des valeurs propres # individus plot(ACM, choix="ind", invisible="var") # individus coloriés selon leurs modalités sur une variable (l'âge) plotellipses(ACM, keepvar="age") # graphique des variables plot(ACM, choix="var") # graphique des modalités plot(ACM, choix="ind", invisible="ind", xlim=c(-1,2), autoLab="yes", cex=0.7) # résultats pour les individus ACM$ind$coord[1:5,] # coordonnées des modalités sur les axes factoriels ACM$var$coord[1:5,] # explication de la modalité par l’axe = carré du coefficient de corrélation de la modalité et de l'axe ACM$var$cos2[1:5,] rowSums(ACM$var$cos2) # somme des cos2 pour l'ensemble des axes = 1 pour chaque modalité sum(colSums(ACM$var$cos2)) # somme des cos2 sum(modal) # contribution de la modalité aux axes factoriels = proportion de l’inertie de la modalité dans l’inertie de l’axe ACM$var$contrib[1:5,] colSums(ACM$var$contrib) # somme des contributions pour chaque colonne = 100 pour chaque axe # On s’intéressera généralement aux modalités dont la contribution est supérieure au poids, # c’est-à-dire dont la coordonnée est supérieure à la racine carrée de l'inertie de l'axe # AFC avec autres packages library(MASS) AFC <- corresp(~ age + rayon, data=afc) print(AFC) plot(AFC) install.packages('ade4') library(ade4) AFC <- dudi.coa(as.data.frame(table(afc$age, afc$rayon)), scannf = FALSE, nf = 2) AFC # --------------------------------------------------------------------------------------------------------- # Chapitre 7 : Analyse factorielle multiple # --------------------------------------------------------------------------------------------------------- library(FactoMineR) data(wine) res <- MFA(wine, group=c(2,5,3,10,9,2), type=c("n",rep("s",5)), ncp=5, name.group=c("orig","olf","vis","olfag","gust","ens"), num.group.sup=c(1,6)) # résultats summary(res) # valeurs propres barplot(res$eig[,1],main="Eigenvalues",names.arg=1:nrow(res$eig)) # explication des axes dimdesc(res) plot(res, choix="group", palette=palette()) plot(res, choix="var", invisible="quanti.sup", hab="group", palette=palette()) plot(res, choix="ind", partial="all", habillage="group", palette=palette()) plot(res, choix="axes", habillage="group", palette=palette()) # --------------------------------------------------------------------------------------------------------- # Chapitre 8 : Classification DBSCAN # --------------------------------------------------------------------------------------------------------- install.packages("factoextra") install.packages("fpc") install.packages("dbscan") install.packages("ggplot2") library(factoextra) data("multishapes", package = "factoextra") table(multishapes[,3]) df <- multishapes[, 1:2] plot(df) library(ggplot2) qplot(df$x, df$y, shape=as.factor(multishapes[,3])) # kmeans set.seed(123) # pour partir des mêmes centres initiaux k <- kmeans(df, centers=6, nstart=10) # 'nstart' feature in kmeans gives the best results out of n trials # (minimizing the total within cluster sum of squares) # taille des classes k$size # croisement entre classes détectées et réelles table(multishapes[,3], k$cluster) # numéro de classe de chaque observation head(k$cluster,10) # centres des classes k$centers # affichage fviz_cluster(k, df, stand = FALSE, ellipse = FALSE, geom = "point") qplot(df$x, df$y, shape=as.factor(k$cluster)) # PAM library(cluster) k <- pam(df, 6) summary(k) plot(k) # taille des classes k$clusinfo # croisement entre classes détectées et réelles table(multishapes[,3], k$cluster) # numéro de classe de chaque observation head(k$cluster,10) # affichage fviz_cluster(k, df, stand = FALSE, ellipse = FALSE, geom = "point") qplot(df$x, df$y, shape=as.factor(k$cluster)) # CLARA library(cluster) k <- clara(df, 6, samples=50) summary(k) plot(k) # taille des classes k$clusinfo # croisement entre classes détectées et réelles table(multishapes[,3], k$cluster) # numéro de classe de chaque observation head(k$cluster,10) # affichage fviz_cluster(k, df, stand = FALSE, ellipse = FALSE, geom = "point") qplot(df$x, df$y, shape=as.factor(k$cluster)) # DBSCAN avec dbscan library(dbscan) dbscan::kNNdistplot(df, k = 5) abline(h = 0.15, lty = 2) k <- dbscan(df, eps = .15, minPts = 5) #k <- dbscan(df, eps = .15, minPts = 5, borderPoints = FALSE) # description de la classification print(k) # nombre des classes table(k$cluster) # croisement entre classes détectées et réelles table(multishapes[,3], k$cluster) # affichage fviz_cluster(k, df, stand = FALSE, ellipse = FALSE, geom = "point") qplot(df$x, df$y, shape=as.factor(k$cluster)) # OPTICS avec dbscan library(dbscan) k <- optics(df, eps = 10, minPts = 10) k # graphique-accessibilité d’OPTICS plot(k) abline(h = 0.2, lty = 2) # extraction du clustering au seuil retenu de distance-accessibilité (eps_cl) res <- extractDBSCAN(k, eps_cl = .2) res plot(res) # noir = bruit hullplot(df, res) # extraction du clustering hiérarchique par la méthode Xi de Ankerst et al. res <- extractXi(k, xi = 0.2) res plot(res) hullplot(df, res) # structure des clusters Xi res$clusters_xi table(res$clusters) table(res$cluster) # croisement entre classes détectées et réelles table(multishapes[,3], res$cluster) # DBSCAN avec fpc library(fpc) set.seed(123) db <- fpc::dbscan(df, eps = 0.15, MinPts = 5) # Plot DBSCAN results plot(db, df, main = "DBSCAN", frame = FALSE) fviz_cluster(db, df, stand = FALSE, ellipse = FALSE, geom = "point") # Print DBSCAN print(db) dbscan::kNNdistplot(df, k = 5) abline(h = 0.15, lty = 2) # --------------------------------------------------------------------------------------------------------- # Chapitre 8 : Mélange de gaussiennes # --------------------------------------------------------------------------------------------------------- # Avec R #----------------------------------------- install.packages("mclust") library(mclust) # Exemple de données : un mélange de deux distributions gaussiennes set.seed(123) data1 <- matrix(rnorm(100, mean = 0, sd = 1), ncol = 2) data2 <- matrix(rnorm(100, mean = 5, sd = 1), ncol = 2) data <- rbind(data1, data2) # Appliquer le modèle de mélange gaussien gmm_model <- Mclust(data) # Résumé du modèle summary(gmm_model) # Visualisation des clusters plot(gmm_model, what = "classification") # Visualisation des probabilités d'appartenance plot(gmm_model, what = "uncertainty") # Analyse de la qualité du clustering # étiquettes de cluster labels <- gmm_model$classification # calcul de la matrice de distances (euclidienne par défaut) dist_matrix <- dist(data) # calcul de l'indice Silhouette sil <- silhouette(labels, dist_matrix) cat("Indice moyen de silhouette :", mean(sil[, 3]), "\n") # Avec PYTHON #----------------------------------------- import numpy as np import matplotlib.pyplot as plt from sklearn.mixture import GaussianMixture # Exemple de données : un mélange de deux distributions gaussiennes np.random.seed(123) data1 = np.random.normal(loc=0, scale=1, size=(100, 2)) data2 = np.random.normal(loc=5, scale=1, size=(100, 2)) data = np.vstack((data1, data2)) # Modèle de mélange gaussien gmm = GaussianMixture(n_components=2, random_state=123) gmm.fit(data) # Prédictions : clusters et probabilités d'appartenance labels = gmm.predict(data) probs = gmm.predict_proba(data) # Visualisation des clusters plt.scatter(data[:, 0], data[:, 1], c=labels, cmap='viridis', s=50, alpha=0.7) plt.scatter(gmm.means_[:, 0], gmm.means_[:, 1], c='red', s=100, marker='x') # Centroides plt.title("Clustering par mélange de gaussiennes (GMM)") plt.xlabel("Feature 1") plt.ylabel("Feature 2") plt.show() # Visualisation des probabilités d'appartenance pour le premier cluster plt.scatter(data[:, 0], data[:, 1], c=probs[:, 0], cmap='coolwarm', s=50, alpha=0.7) plt.title("Probabilités d'appartenance au premier cluster") plt.xlabel("Feature 1") plt.ylabel("Feature 2") plt.colorbar(label="Probabilité") plt.show() # Analyse de la qualité du clustering from sklearn.metrics import silhouette_score silhouette = silhouette_score(data, labels) print(f"Indice Silhouette : {silhouette:.3f}") # Distribution des points par cluster unique, counts = np.unique(labels, return_counts=True) print("\nRépartition des points par cluster :") for u, c in zip(unique, counts): print(f" Cluster {u} : {c} points ({100*c/len(data):.1f}%)") # --------------------------------------------------------------------------------------------------------- # Chapitre 8 : Classification par agrégation de similarités # --------------------------------------------------------------------------------------------------------- install.packages('amap') library(amap) credit <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data", header = FALSE) for (i in 1:21) credit[,i] <- factor(credit[,i]) str(credit) creditn <- matrix(c(lapply(credit,as.integer),recursive=T), ncol=21) head(creditn) str(creditn) matrice <- diss(creditn) matrice[1:10,1:10] pop(matrice) # calcul très long ! etude <- read.csv2("C:/Users/.../enquete.csv") str(etude) head(etude) summary(etude) etuden <- matrix(c(lapply(etude,as.integer),recursive=T), ncol=6) matrice <- diss(etuden) dim(matrice) system.time(pop(matrice)) # --------------------------------------------------------------------------------------------------------- # Chapitre 8 : Classification de variables # --------------------------------------------------------------------------------------------------------- etude <- read.csv2("C:/Users/.../etude_de_cas.csv", stringsAsFactors = TRUE) str(etude) head(etude) install.packages('ClustOfVar') library(ClustOfVar) # classification de quelques variables quantitatives varquanti <- c("age", "anciennete", "relation", "nbpoints", "nbproduits", "nbachats", "revenus", "abonnement1", "evolconsom", "utilcredit") tree <- hclustvar(etude[names(etude) %in% varquanti]) # dendrogramme plot(tree) # stabilité stab <- stability(tree, B=100) plot(stab, main="Stabilité des partitions") plot(stab, nmax=7) # coupure du dendrogramme part_hier <- cutreevar(tree, 3, matsim = TRUE) print(part_hier) summary(part_hier) part_hier$var # coeff de corrélation entre chaque variable quali et chaque variable synthétique d'un cluster cor(part_hier$scores) head(part_hier$scores, 10) part_hier$cluster part_hier$size part_hier$sim # corrélations des variables à l'intérieur de chaque classe # homogénéité part_hier$wss # gain d'homogénéité sum(part_hier$wss) 100*(sum(part_hier$wss) - cutreevar(tree, 1)$wss) / (sum(part_hier$size) - cutreevar(tree, 1)$wss) part_hier$E 100*sum(part_hier$wss) / sum(part_hier$size) # classification des variables quantitatives et qualitatives # séparation des variables qualis et quantis classes <- as.character(sapply(etude, class)) (varquali <- which(classes=="factor")) (varquanti <- which(classes!="factor")) (varquanti <- which(classes!="factor")[-1]) # suppression de l'identifiant client (variable numérique) # on appelle la fonction hclustvar de CAH avec les variables # - quantis en 1er argument # - qualis en 2e argument tree <- hclustvar(etude[,varquanti], etude[,varquali]) # gain d'homogénéité de la classification en fonction du nb de classes de variables gain <- sapply(2:6, function(n) (cutreevar(tree,n)$E)) barplot(gain) # on appelle la fonction kmeansvar de kmeans avec les variables # - quantis en 1er argument # - qualis en 2e argument set.seed(123) #part_km <- kmeansvar(etude[names(etude) %in% varquanti], init=4, nstart=5) part_km <- kmeansvar(etude[,varquanti], etude[,varquali], init=3, nstart=10) part_km$var # homogénéité des classes = somme des carrés des coefficients de corrélation # des variables de la classe avec la composante principale de la classe part_km$E part_hier$E # gain d'homogénéité de la classification en fonction du nb de classes de variables gain <- sapply(2:6, function(x) {set.seed(123) ; kmeansvar(etude[,varquanti], etude[,varquali], init=x, nstart=10)$E}) # --------------------------------------------------------------------------------------------------------- # Chapitre 10 : Régression linéaire clusterwise # --------------------------------------------------------------------------------------------------------- set.seed(2) x1 <- rnorm(100) set.seed(3) y1 <- x1 + rnorm(100, sd=0.5) set.seed(5) y2 <- - x1 + rnorm(100, sd=0.5) x <- c(x1,x1) y <- c(y1,y2) modele <- lm(y ~ x) summary(modele) plot(x,y) abline(modele, lty=3) install.packages('flexmix') library(flexmix) clw <- flexmix(y ~ x, k=2) summary(clw) parameters(clw, component = 1) parameters(clw, component = 2) summary(refit(clw)) plot(clw) plot(x, y, pch=clusters(clw), col=clusters(clw)) plot(x, y, pch=clusters(clw), col=c("black", "grey")[clusters(clw)]) abline(parameters(clw, component = 1)[1:2], lty=3, col="black") abline(parameters(clw, component = 2)[1:2], lty=4, col="grey") # --------------------------------------------------------------------------------------------------------- # Chapitre 12 : Régression linéaire # --------------------------------------------------------------------------------------------------------- # 1er jeu de données model <- lm(dist ~ speed, data=cars) # 2e jeu de données data(anscombe) model <- lm(y1 ~ x1+x4, data=anscombe) # modèle summary(model) predict(model, newdata = data.frame(speed=7, dist=2), interval='prediction', level = 0.95) # intervalles de confiance des prédictions et de la moyenne pred.int <- predict(model, interval = "prediction") df <- cbind(cars, pred.int) library(ggplot2) p <- ggplot(df, aes(speed, dist)) + geom_point() + stat_smooth(method = lm) p + geom_line(aes(y = lwr), color="red", linetype="dashed") + geom_line(aes(y = upr), color="red", linetype="dashed") # exemple de la consommation de chauffage chauffage <- data.frame(consommat=c(1042,1377,622,154,357,874,1388,1138,900,460,119,770,1670,1223,199), temperat=c(4.4,-2.8,4.4,22.8,17.8,1.1,-12.8,-13.3,-5.0,17.2,18.3,5.0,-6.1,3.3,14.4), isolatio=c(7.6,7.6,25.4,15.2,15.2,15.2,15.2,25.4,25.4,7.6,25.4,15.2,7.6,7.6,25.4)) library(car) scatter3d(chauffage$consommat, chauffage$temperat, chauffage$isolatio, fit="linear", residuals=TRUE, bg="white", axis.scales=FALSE, grid=TRUE, xlab="consommat",ylab="temperat", zlab="isolatio") chauffage.lm <- lm(consommat~temperat+isolatio, data=chauffage) summary(chauffage.lm) par(mfrow=c(2,2)) plot(chauffage.lm) # test de Shapiro (on veut une p-value > 5% pour ne pas rejeter l'hypothèse H0 selon laquelle les résidus sont issus d'une population normalement distribuée) shapiro.test(model$residuals) # test de Breusch-Pagan (on veut une p-value > 5% pour ne pas rejeter l'hypothèse H0 d’homoscédasticité des résidus) library(lmtest) bptest(model) install.packages('skedastic') library(skedastic) breusch_pagan(model) # à la main bpmodel <- lm(residuals(model)^2 ~ x1+x4, data=anscombe) (studentized_bp <- nrow(anscombe) * summary(bpmodel)$r.squared) pchisq(studentized_bp, df=2, lower.tail = FALSE) # test de White library(skedastic) # avec interactions (5 degrés de liberté) white(chauffage.lm, interactions = TRUE) # sans interactions (4 degrés de liberté) white(chauffage.lm, interactions = FALSE) # à la main rescarre <- residuals(chauffage.lm)^2 # avec interactions (5 degrés de liberté) white <- lm(rescarre ~ (temperat*isolatio) + poly(temperat,2) + poly(isolatio,2), data = chauffage) white.stat <- summary(white)$r.squared * nrow(chauffage) c(white.stat, pchisq(white.stat, df = 5, lower = FALSE)) # sans interactions (4 degrés de liberté) white <- lm(rescarre ~ poly(temperat,2) + poly(isolatio,2), data = chauffage) white.stat <- summary(white)$r.squared * nrow(chauffage) c(white.stat, pchisq(white.stat, df = 4, lower = FALSE)) # test de Durbin-Watson (on veut une p-value > 5% pour ne pas rejeter l'hypothèse H0 d'absence d'autocorrélation) library(lmtest) dwtest(model) # à la main residus <- residuals(model) diffs <- diff(residus) (dw <- sum(diffs^2)/sum(residus^2)) # Cochrane-Orcutt moncochrane <- function(m) { co.rho <- coefficients(lm(m$residuals[-length(m$residuals)] ~ 0 + m$residuals[-1])) yt <- y - co.rho*lag(y) xt <- x - co.rho*lag(x) return(coefficients(lm(yt ~ xt))/c((1-co.rho),1)) } moncochrane(model) # colinéarité colin <- data.frame(X1=c(1,2,3,4), X2=c(1.01,1.99,3.01,3.99), Y=c(16,34,44,46), Z=c(17,34,44,46)) cor(colin) colin1.lm <- lm(Y~X1+X2, data=colin) colin2.lm <- lm(Z~X1+X2, data=colin) summary(colin1.lm) summary(colin2.lm) library(car) scatter3d(colin$X1, colin$Y, colin$X2, fit="linear", residuals=TRUE,bg="white", axis.scales=FALSE, grid=TRUE, ellipsoid=FALSE, xlab="X1",ylab="Y", zlab="X2") scatter3d(colin$X1, colin$Z, colin$X2, fit="linear", residuals=TRUE,bg="white", axis.scales=FALSE, grid=TRUE, ellipsoid=FALSE, xlab="X1", ylab="Z", zlab="X2") # --------------------------------------------------------------------------------------------------------- # Chapitre 12 : Régression PLS # --------------------------------------------------------------------------------------------------------- library(sas7bdat) chenilles <- read.sas7bdat("C:/Users/.../chenilles.sas7bdat") print(chenilles) lin <- lm(log~X1+X2+X4+X5, data=chenilles) summary(lin) library(MASS) ridge <- lm.ridge(log~X1+X2+X4+X5, data=chenilles, lambda=seq(0,1,0.05)) plot(ridge) install.packages('pls') library(pls) simpls <- mvr(log~X1+X2+X4+X5, data=chenilles, ncomp=3, method="simpls", scale=TRUE) simpls <- plsr(log~X1+X2+X4+X5, data=chenilles, validation="CV") simpls <- plsr(log~X1+X2+X4+X5, data=chenilles, ncomp=4, scale=F, validation="CV") b <- coef(simpls, ncomp=1, intercept=FALSE) summary(simpls) coef(simpls) loadings(simpls) plot(simpls) plot(RMSEP(simpls)) sapply(1:4, function(i) coef(simpls, ncomp=i, intercept=TRUE)) sapply(1:4, function(i) coef(mvr(log~X1+X2+X4+X5, data=chenilles, ncomp=i, method="simpls", scale=FALSE), intercept=TRUE)) t(sapply(1:4, function(i) coef(plsr(log~X1+X2+X4+X5, data=chenilles, ncomp=i, method="simpls", scale=TRUE), intercept=TRUE))) # par défaut, plsr centre les variables mais ne les réduit pas ; il faut utiliser l'argument SCALE=TRUE validation="CV" msepcv <- MSEP(simpls, estimate=c("train","CV")) plot(msepcv, lty=1, type="l", legendpos="topright", main="") which.min(msepcv$val["CV",,])-1 # régression PLS avec package plsdepot d'après le livre de Michel Tenenhaus : "La Régression PLS: Théorie et Pratique" install.packages('plsdepot') library(plsdepot) names(chenilles) # méthode PLS1 (Partial Least Squares Regression for the univariate case) du package plsdepot pls1 <- plsreg1(chenilles[, c(1,2,4,5)], chenilles$log, comps = 4) pls1 <- plsreg1(scale(chenilles[, c(1,2,4,5)]), scale(chenilles$log), comps = 4) plot(pls1) pls1$x.scores pls1$std.coefs pls1$reg.coefs pls1$R2 pls1$R2Xy pls1$Q2 t(sapply(2:4, function(i) plsreg1(chenilles[, c(1,2,4,5)], chenilles$log, comps = i)$reg.coefs)) t(sapply(2:4, function(i) plsreg1(chenilles[, c(1,2,4,5)], chenilles$log, comps = i)$std.coefs)) # graphique des coefficients en fonction du nombre de facteurs res <- data.frame(sapply(2:4, function(i) plsreg1(chenilles[, c(1,2,4,5)], chenilles$log, comps = i)$reg.coefs[-1])) colnames(res) <- c("2", "3", "4") res$coef <- rownames(res) install.packages('reshape') library(reshape) dm <- melt(res, id="coef") colnames(dm) <- c("coef", "nbfact", "coefficients") library(ggplot2) ggplot(dm, aes(x=nbfact, y=coefficients, linetype=coef, group=coef)) + geom_line() # méthode NIPALS (Non-linear Iterative Partial Least Squares) du package plsdepot ni <- nipals(chenilles, comps = 4) ni$values ni$scores ni$loadings # --------------------------------------------------------------------------------------------------------- # Chapitre 13 : Régression logistique # --------------------------------------------------------------------------------------------------------- # Figure 13.1 : Comparaison des régressions linéaire et logistique set.seed(123) n <- 1000 sd <- 1 x <- c(rnorm(n, sd=sd), 1+rnorm(n,sd=sd)) y <- c(rep(0,n), rep(1,n)) plot(y~x) # prédiction brutale m1 <- mean(x[y==0]) m2 <- mean(x[y==1]) m <- mean(c(m1,m2)) if(m1m2) a <- 1 if(m1==m2) a <- .5 lines( c(min(x),m,m,max(x)), c(a,a,1-a,1-a), lty=1, lwd=1) # régression linéaire abline(lm(y~x), lty=4, lwd=1) # régression logistique logit xp <- seq(min(x), max(x), length=200) mod <- glm(y~x, family=binomial(link = "logit")) summary(mod) yp <- predict(mod, data.frame(x=xp), type='response') lines(xp, yp, lty=5, lwd=2, col=gray(0.5)) # régression logistique probit xp <- seq(min(x), max(x), length=200) mod <- glm(y~x, family=binomial(link = "probit")) summary(mod) yp <- predict(mod, data.frame(x=xp), type='response') lines(xp, yp, lty=6, lwd=1, col=gray(0.8)) # courbe théorique curve(dnorm(x,1,1)*.5/(dnorm(x,1,1)*.5+dnorm(x,0,1)*(1-.5)), add=T, lty=3, lwd=2, col="black") legend( .95*par('usr')[1]+.05*par('usr')[2], .9, #par('usr')[4], c('prédiction simpliste', "régression linéaire", "régression logit", "régression probit", "courbe théorique"), lty=c(1,4,5,6,3), lwd=c(1,1,2,1,2), col=c("black","black",gray(0.5),gray(0.8),"black")) title(main="Comparaison des régressions linéaire et logistique") # Exemple CHD (Coronary Heart Disease)de Hosmer et Lemeshow # https://rdrr.io/cran/aplore3/man/chdage.html install.packages('aplore3') library(aplore3) head(chdage,10) summary(chdage) # figure 1.1 p. 5 plot(as.integer(chd)-1 ~ age, pch = 20, main = "Figure 1.1 p. 5", ylab = "Coronary heart disease", xlab = "Age (years)", data = chdage) # table 1.2 with(chdage, addmargins(table(agegrp))) with(chdage, addmargins(table(agegrp, chd))) (Means <- with(chdage, tapply(as.integer(chd)-1, list(agegrp), mean))) # figure 1.2 p. 6 midPoints <- c(24.5, seq(32, 57, 5), 64.5) plot(midPoints, Means, pch = 20, ylab = "Coronary heart disease (mean)", xlab = "Age (years)", ylim = 0:1, main = "Figure 1.2 p. 6") lines(midPoints, Means) # table 1.3 summary( mod1.3 <- glm( chd ~ age, family = binomial, data = chdage )) # table 1.4 vcov(mod1.3) # computing OddsRatio and confidence intervals for age ... exp(coef(mod1.3))[-1] exp(confint(mod1.3))[-1, ] # Comparaison entre méthode de Newton et méthode de scoring de Fisher # jeu de données set.seed(123) n <- 10000 p <- 5 A <- matrix (runif(n*p), n, p) A <- data.frame(A, y = sample(c(0,1), n, replace = T)) Y <- A[,"y"] X <- as.matrix(cbind(1, A[,!names(A) %in% "y"])) # régression logistique avec glm system.time(L <- glm(y ~., data=A, family=binomial(link="logit"))) summary(L) # régression logistique avec la méthode de Newton "à la main" maxiter <- 10 beta <- lm(Y~0+X)$coefficients for (s in 1:maxiter){ pi <- 1/(1+exp(-X%*%beta)) gradient <- t(X)%*%(Y-pi) omega <- matrix(0,nrow(X),nrow(X)) diag(omega) <- (pi*(1-pi)) hessian <- -t(X)%*%omega%*%X if (sum(abs(solve(hessian)%*%gradient)) < (p*0.000001)) break beta <- beta - solve(hessian)%*%gradient } sd <- sqrt(diag(solve(-hessian))) # p-value associée au test de Student modele <- data.frame(beta, sd, beta/sd, 2*(1-pnorm(abs(beta/sd)))) colnames(modele) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)") modele cat("\n", "Convergence en ", s-1 ," itérations","\n") # régression logistique IRLS (méthode de Fisher) IRLS <- function(x,y) { change <- Inf beta.old <- lm(y~0+x)$coef while(change > 1e-8) { pi <- 1/(1+exp(-x%*%beta.old)) weight <- pi*(1-pi) z <- ((y-pi)/weight) + log(pi/(1-pi)) beta.new <- lm(z~0+x, weights=weight)$coef change <- sqrt(sum((beta.new - beta.old)^2)) beta.old <- beta.new } s <- coef(summary(lm(z~0+x, weights=weight)))[,"Std. Error"] sd <- s/sqrt(sum(((y-pi)^2)/weight)/(n-p-1)) # p-value associée au test de Student modele <- data.frame(beta.new, sd, beta.new/sd, 2*(1-pnorm(abs(beta.new/sd)))) colnames(modele) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)") print(modele) } system.time(IRLS(X,Y)) # régression logistique IRLS (sans le calcul des erreurs-types des coefficients) IRLS <- function(x,y) { change <- Inf beta.old <- lm(y~0+x)$coef while(change > 1e-8) { pi <- 1/(1+exp(-x%*%beta.old)) weight <- pi*(1-pi) z <- ((y-pi)/weight) + log(pi/(1-pi)) beta.new <- lm(z~0+x, weights=weight)$coef change <- sqrt(sum((beta.new - beta.old)^2)) beta.old <- beta.new } #print(beta.new) } system.time(IRLS(X,Y)) # régression logistique IRLS (lm() remplacée par calcul matriciel direct) library(car) IRLS0 <- function(x,y) { change <- Inf beta.old <- solve(crossprod(x), crossprod(x,y)) while(change > 1e-8) { pi <- 1/(1+exp(-x%*%beta.old)) weight <- pi*(1-pi) z <- ((y-pi)/weight) + log(pi/(1-pi)) beta.new <- solve(wcrossprod(x,w=as.numeric(weight)), wcrossprod(x,z,w=as.numeric(weight))) change <- sqrt(sum((beta.new - beta.old)^2)) beta.old <- beta.new } #print(beta.new) } system.time(IRLS0(X,Y)) # méthode de Newton Newton <- function(x,y) { change <- Inf beta.old <- lm(y~0+x)$coef while(change > 1e-8) { pi <- 1/(1+exp(-x%*%beta)) gradient <- t(x)%*%(y-pi) omega <- matrix(0,nrow(x),nrow(x)) diag(omega) <- (pi*(1-pi)) hessian <- -t(x)%*%omega%*%x change <- sqrt(sum((solve(hessian)%*%gradient)^2)) beta <- beta - solve(hessian)%*%gradient } #print(beta) } system.time(Newton(X,Y)) install.packages('microbenchmark') library(microbenchmark) microbenchmark(Newton=Newton(X,Y), IRLS=IRLS(X,Y), IRLS0=IRLS0(X,Y), GLM=glm(y ~., data=A, family=binomial(link="logit")), times=100, unit="ms") # Tests statistiques de la régression logistique # R2 de Cox-Snell (r2cs <- 1- (exp(-summary(logit)$null.deviance/(1+logit$df.null))/exp(-summary(logit)$deviance/(1+logit$df.null)))) # R2 de Nagelkerke (r2max <- 1- (exp(-summary(logit)$null.deviance/(1+logit$df.null)))) (r2nag <- r2cs/r2max) # R2 de McFadden (r2mcfd <- 1- summary(logit)$deviance/summary(logit)$null.deviance) logLik(logit) # Test de Hosmer-Lemeshow install.packages('aplore3') library(aplore3) model <- glm(chd~age, data=chdage, family=binomial(link = "logit")) summary(model) # implémentation 1 du test de Hosmer et Lemeshow install.packages('glmtoolbox') library(glmtoolbox) hltest(model) # implémentation 2 du test de Hosmer et Lemeshow install.packages('generalhoslem') library(generalhoslem) (hl <- logitgof(chdage$chd, fitted(model))) hl$observed sum(hl$observed) hl$expected sum(hl$expected) # découpage de la fonction hldf <- data.frame(cbind(hl$observed,hl$expected)) hldf # calcul du chi2 (hlt <- sum((hldf$y1-hldf$yhat1)^2/hldf$yhat1)+sum((hldf$y0-hldf$yhat0)^2/hldf$yhat0)) pchisq(hlt, df = 8, lower.tail = FALSE) # découpage alternatif (= version 2) hldf <- data.frame(y0=c(9,9,8,8,7,5,5,3,2,1),y1=c(1,1,2,3,4,5,5,10,8,4), yhat0=c(9.213,8.657,8.095,8.037,6.947,5.322,4.2,3.736,2.134,0.661),yhat1=c(0.787,1.343,1.905,2.963,4.053,4.678,5.8,9.264,7.866,4.339)) hldf (hlt <- sum((hldf$y1-hldf$yhat1)^2/hldf$yhat1)+sum((hldf$y0-hldf$yhat0)^2/hldf$yhat0)) pchisq(hlt, df = 8, lower.tail = FALSE) # implémentation 3 du test de Hosmer et Lemeshow (Modified Hosmer-Lemeshow Test for Large Samples) # https://rdrr.io/github/gnattino/largesamplehl/man/hltest.html install.packages('devtools') library(devtools) install_github("gnattino/largesamplehl") library(largesamplehl) hltest(model) # valeur de epsilon0 sqrt((qchisq(0.95, df = 8)-8)/1e6) sqrt((qchisq(0.05, df = 8, lower.tail = FALSE)-8)/1e6) # sensibilité du test de Hosmer et Lemeshow à l'effectif chdage_rep <- do.call(rbind, replicate(10, chdage, simplify = FALSE)) model <- glm(chd~age, data=chdage_rep, family=binomial(link = "logit")) summary(model) hltest(model) logitgof(chdage_rep$chd, fitted(model)) logitgof(chdage_rep$chd, fitted(model), g=15) # --------------------------------------------------------------------------------------------------------- # Chapitre 14 : Séparation complète en régression logistique # --------------------------------------------------------------------------------------------------------- X <- c(9,10,11,11,12,14,14,15,16,18,12,13,15,17,19,18,20,21,22,23) Y <- c(42,63,54,68,48,37,55,60,52,62,15,9,20,5,19,44,34,11,40,23) C1 <- c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1) C2 <- c(0,0,0,0,0,0,0,0,1,0,1,1,1,1,1,1,1,1,1,1) C3 <- c(0,0,0,1,0,0,0,0,1,0,1,1,1,1,1,1,1,1,1,1) plot(X, Y, col=factor(C1), cex=2, pch=16) modele1 <- glm(C1~X+Y, family=binomial(link = "logit"), control = list(maxit = 20)) summary(modele1) # séparation incomplète plot(X, Y, col=factor(C2), cex=2, pch=16) modele2 <- glm(C2~X+Y, family=binomial(link = "logit")) summary(modele2) # séparation incomplète plot(X, Y, col=factor(C3), cex=2, pch=16) modele3 <- glm(C3~X+Y, family=binomial(link = "logit")) summary(modele3) # visualisation des prédictions pred_logistic = function (x, y) { return(predict(modele1, newdata=data.frame(X=x, Y=y), type="response")) } pred_logistic2 = function (x, y) { return(predict(modele2, newdata=data.frame(X=x, Y=y), type="response")) } vx <- seq (min(X)-1,max(X)+1, length=50) vy <- seq (min(Y)-1,max(Y)+1, length=50) par(mfrow=c(1,2)) image(vx, vy , outer(vx, vy, pred_logistic), col=gray((1:32)/32), xlab="X", ylab="Y") points(X, Y, col=c('red', 'blue')[as.numeric(C1)+1], cex=2, pch=16) mtext("Séparation complète", side=1, line=4, cex=1) image(vx, vy , outer(vx, vy, pred_logistic2), col=gray((1:32)/32), xlab="X", ylab="Y") points(X, Y, col=c('red', 'blue')[as.numeric(C1)+1], cex=2, pch=16) mtext("Séparation incomplète", side=1, line=4, cex=1) par(mfrow=c(1,1)) # régression logistique de Firth (1993), "Bias reduction of maximum likelihood estimates", Biometrika, 80,1 install.packages('logistf') library(logistf) firth <- logistf(C1~X+Y, family=binomial(link = "logit")) summary(firth) head(predict(firth, type="response")) # approche bayésienne utilisant l'information a priori via l'algorithme EM install.packages('arm') library(arm) bayes <- bayesglm(C1 ~ X + Y, prior.scale=Inf, prior.df=Inf, family="binomial") # logit classique bayes <- bayesglm(C1 ~ X + Y, prior.scale=2.5, prior.df=1, family="binomial") # prior Cauchy avec e-t 2,5 bayes <- bayesglm(C1 ~ X + Y, prior.scale=2.5, prior.df=Inf, family="binomial") # prior Normal avec e-t 2,5 display(bayes) summary(bayes) head(predict(bayes, type="response")) # visualisation des prédictions pred_logistic = function (x, y) { return(predict(modele1, newdata=data.frame(X=x, Y=y), type="response")) } pred_firth = function (x, y) { return(predict(firth, newdata=data.frame(X=x, Y=y), type="response")) } pred_bayes = function (x, y) { return(predict(bayes, newdata=data.frame(X=x, Y=y), type="response")) } vx <- seq (min(X)-1,max(X)+1, length=50) vy <- seq (min(Y)-1,max(Y)+1, length=50) par(mfrow=c(1,3)) image(vx, vy , outer(vx, vy, pred_logistic), col=gray((1:32)/32), xlab="X", ylab="Y") points(X, Y, col=c('red', 'blue')[as.numeric(C1)+1], cex=2, pch=16) mtext("Régression logistique", side=1, line=4, cex=1) image(vx, vy , outer(vx, vy, pred_firth), col=gray((1:32)/32), xlab="X", ylab="Y") points(X, Y, col=c('red', 'blue')[as.numeric(C1)+1], cex=2, pch=16) mtext("Régression de Firth", side=1, line=4, cex=1) image(vx, vy , outer(vx, vy, pred_bayes), col=gray((1:32)/32), xlab="X", ylab="Y") points(X, Y, col=c('red', 'blue')[as.numeric(C1)+1], cex=2, pch=16) mtext("Régression bayésienne", side=1, line=4, cex=1) par(mfrow=c(1,1)) # section 13.1.14 : régression logistique sur jeu de données Titanic install.packages('explore') library(explore) titanic <- as.data.frame(use_data_titanic(count = FALSE)) # remplacement des chaînes de caractère par des numéros de classe titanic[] <- lapply(titanic, function(x) as.integer(as.factor(x))-1) # on veut que la classe "crew" ait le numéro 0, la 1ère classe le numéro 1, etc. titanic$Class <- (1+titanic$Class)%%4 titanic$Age <- 1-titanic$Age head(titanic) attach(titanic) table (Class,Survived) prop.table(table(Class,Survived),1) # régression logistique sur variables quantitatives summary(glm(Survived~Class+Age+Sex, data=titanic, family=binomial(link = "logit"))) summary(titanic) for (i in 1:4) titanic[,i] <- factor(titanic[,i]) summary(titanic) titanic$Class <- relevel(titanic$Class, ref="3") # régression logistique sur variables qualitatives summary(glm(Survived~Class+Age+Sex, data=titanic, family=binomial(link = "logit"))) # --------------------------------------------------------------------------------------------------------- # Chapitre 13 : Régression logistique pénalisée # --------------------------------------------------------------------------------------------------------- # régression elasticnet (méthode de descente des coordonnées) library(glmnet) # jeu de données credit <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data", header=FALSE, stringsAsFactors=TRUE) # noms des variables colnames(credit) <- 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") credit$Etranger <- NULL # variable à expliquer credit$Cible[credit$Cible == 1] <- 0 # crédits OK credit$Cible[credit$Cible == 2] <- 1 # crédits KO credit$Cible <- factor(credit$Cible) # discrétisation des variables continues credit$Age <- cut(credit$Age,c(0,25,Inf),right=TRUE) #intervalle semi-fermé à droite credit$Duree_credit <- cut(credit$Duree_credit,c(0,15,36,Inf),right=TRUE) credit$Montant_credit <- cut(credit$Montant_credit,c(0,4000,Inf),right=TRUE) # transformation en facteurs des variables discrètes credit[["Taux_effort"]] <- factor(credit[["Taux_effort"]]) credit[["Anciennete_domicile"]] <- factor(credit[["Anciennete_domicile"]]) credit[["Nb_credits"]] <- factor(credit[["Nb_credits"]]) credit[["Nb_pers_charge"]] <- factor(credit[["Nb_pers_charge"]]) # tirage aléatoire stratifié sans remise library(sampling) set.seed(235) id <- strata(credit, stratanames="Cible", size=c(sum(credit$Cible==0)*2/3,sum(credit$Cible==1)*2/3), method="srswor", description=T)$ID_unit # échantillons d'apprentissage et de validation train <- credit[id,] valid <- credit[-id,] # 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"] penfit = cv.glmnet(X, Y, alpha=0, family = "binomial", type="auc", nlambda=100) plot(penfit) # graphique penfit$lambda[1] # plus petit lambda annulant tous les coefficients penfit$lambda[100] # lambda précédent divisé par 10 000 fit = glmnet(X, Y, alpha=0, family = "binomial", lambda=seq(penfit$lambda[100], penfit$lambda[1], length=10000), standardize = FALSE) plot(fit, xvar="lambda", sign.lambda = -1) # --------------------------------------------------------------------------------------------------------- # Chapitre 14 : Méthode MARS (Multivariate Adaptive Regression Splines) # --------------------------------------------------------------------------------------------------------- install.packages('earth') library(earth) data(ozone1) mars <- earth(O3 ~ temp, data = ozone1) plot(mars) plotmo(mars) a <- earth(Volume ~ Girth, data=trees, linpreds=1) a <- earth(Volume ~ Girth, data=trees, linpreds=1, glm=list(family=gaussian)) summary(a, details=TRUE) # details=TRUE pour avoir les p-values des coefficients avec l'option "glm" ci-dessus plot(a) plotmo(a) plotmo(a, caption="", col.response="red", clip=0, xlab="x", ylab="y", main="") plotmo(a, col.response="grey", xlab="x", ylab="y", main="Modèle linéaire") b <- earth(Volume ~ Girth, data=trees) b <- earth(Volume ~ Girth, data=trees, glm=list(family=gaussian)) summary(b, details=TRUE) # details=TRUE pour avoir les p-values des coefficients avec l'option "glm" ci-dessus plot(b) par(mfrow=c(2,2)) plot(b$glm.list[[1]]) summary(b, decomp="none") summary(lm(trees$Volume ~ b$bx - 1)) plotmo(b, caption="", col.response="red", clip=0, xlab="x", ylab="y", main="") plotmo(b, col.response="grey", xlab="x", ylab="y", main="Modèle MARS") # résidus head(b$residuals) shapiro.test(b$residuals) library(lmtest) bptest(b) dwtest(b) # modèle linéaire classique l <- lm(Volume ~ Girth, data=trees) summary(l) # pour comparaison, une autre fonction de régression linéaire par morceaux : # piecewise regression install.packages('segmented') library(segmented) # fit simple linear regression model fit <- lm(Volume ~ Girth, data=trees) # fit piecewise regression model to original model segmented.fit <- segmented(fit, seg.Z = ~Girth) # view summary of segmented model summary(segmented.fit) # plot original data plot(trees$Girth, trees$Volume, col='red') # add segmented regression model plot(segmented.fit, add=TRUE) # --------------------------------------------------------------------------------------------------------- # Chapitre 15 : Figures sur les SVM # --------------------------------------------------------------------------------------------------------- # Figure 15.1 pour livre (Séparation correcte (B2) et séparation optimale (B1)) library(e1071) x1s <- c(.5,1,1,2,3,3.5,1,3.5,4,5,5.5,2) x2s <- c(3.5,1,2.5,2,1,1.2,5.8,3,4,5,4,4.5) ys <- c(rep(+1,6), rep(-1,6)) my.data <- data.frame(x1=x1s, x2=x2s, type=as.factor(ys)) my.data plot(my.data[,-3], pch=ifelse(ys==1,1,15), cex = 1, xlim=c(-1,6), ylim=c(-1,6)) svm.model <- svm(type ~ ., data=my.data, type='C-classification', kernel='linear', scale=FALSE) #points(my.data[svm.model$index,c(1,2)],col="blue",cex=2) # vecteurs supports points(my.data[svm.model$index,c(1,2)], col="grey30", cex=3) # vecteurs supports x1min = min(x1s); x1max = max(x1s); x2min = min(x2s); x2max = max(x2s); coef1 = sum(svm.model$coefs*x1s[svm.model$index]) coef2 = sum(svm.model$coefs*x2s[svm.model$index]) lines(c(x1min,x1max), (svm.model$rho-coef1*c(x1min, x1max))/coef2) lines(c(x1min,x1max), (svm.model$rho+1-coef1*c(x1min, x1max))/coef2, lty=2) lines(c(x1min,x1max), (svm.model$rho-1-coef1*c(x1min, x1max))/coef2, lty=2) lines(c(0,5.5),c(4,2.4),lty=3) lines(c(0,5.5),c(3.7,2.1),lty=3) lines(c(0,5.5),c(3.85,2.25),lty=4) legend("bottomleft", c("B1 : Séparation optimale","B2 : Séparation correcte"), lty=c(2,3), col=c("black","black"), lwd=2, seg.len=4) text(0,4.7,"B1") text(-0.5,4,"B2") # Figure 15.3 illustrant le kernel trick library(lattice) install.packages('gridExtra') library(gridExtra) n <- 100 set.seed(235) x <- runif(n,-1,1) set.seed(456) y <- runif(n,-1,1) #classe <- ifelse(x^2+y^2<0.5,"a","b") classe <- ifelse((x/2)^2+y^2<0.5,"a","b") ellipse.x <- seq(-1.5,1.5,by=0.02) ellipse.y <- sqrt(0.5-((ellipse.x)/2)^2) x1 <- c(x,ellipse.x,ellipse.x) x2 <- c(y,ellipse.y,-ellipse.y) groupe <- c(classe,rep("c",2*length(ellipse.x))) z1 <- x1^2 z2 <- sqrt(2)*x1*x2 z3 <- x2^2 p1 <- xyplot(x2 ~ x1, scales=list(draw=F), groups=groupe, par.settings=list(superpose.symbol=list(pch=c(15,2,8),cex=c(1,1,0.3),col='black'))) p2 <- cloud(z3 ~ z1+z2, groups=groupe, par.settings=list(superpose.symbol=list(pch=c(15,2,8),cex=c(1,1,0.3),col='black'))) grid.arrange(p1, p2, ncol=2) # --------------------------------------------------------------------------------------------------------- # Chapitre 15 : Reconnaissance de chiffres manuscrits avec des SVM (base MNIST) # --------------------------------------------------------------------------------------------------------- # les images sont de taille 28 x 28 pixels # le fichier train.csv a 785 colonnes : # - label = chiffre # - pixel0 à pixel783 sont les 784 = 28 x 28 pixels de l'image du chiffre # les valeurs des pixels sont des niveaux de gris entre 0 et 255 # réf : http://yann.lecun.com/exdb/mnist/ # données directement disponibles au format csv : http://www.pjreddie.com/projects/mnist-in-csv/ train <- read.csv("C:/Users/.../mnist_train.csv", header=TRUE) dim(train) test <- read.csv("C:/Users/.../mnist_test.csv", header=TRUE) dim(test) # Modélisation sur facteurs ACP avec sortie softmax set.seed(235) trainDataset <- train[sample(nrow(train)),] # mélange des données X <- trainDataset[,-1] Y <- factor(trainDataset[,1]) Xreduced <- X/255 # ACP Xcov <- cov(Xreduced) pcaX <- prcomp(Xcov) # matrices d'apprentissage en entrée du réseau de neurones nbaxes <- 58 Xfinal <- as.matrix(Xreduced) %*% pcaX$rotation[,1:nbaxes] dim(Xfinal) rm(Xreduced) gc() # données de test transformées par ACP testreduced <- test[,-1]/255 testreduced <- as.matrix(testreduced) %*% pcaX$rotation[,1:nbaxes] # SVM sur axes factoriels library(e1071) # SVM avec noyau linéaire system.time(svmlin <- svm(Xfinal, Y, kernel="linear", probability=TRUE, cost=0.1)) summary(svmlin) svmlin$rho pred.svm <- predict(svmlin, newdata=testreduced, decision.values=T) # taux d'erreur en test sum(pred.svm != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.svm) # SVM avec noyau radial system.time(svmrad <- svm(Xfinal, Y, kernel="radial", probability=TRUE, gamma=0.01, cost=10)) summary(svmrad) svmrad$rho pred.svm <- predict(svmrad, newdata=testreduced, decision.values=TRUE) # taux d'erreur en test sum(pred.svm != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.svm) # SVM avec noyau polynomial system.time(svmpol <- svm(Xfinal, Y, kernel="polynomial", degree=2, coef0=0.2, gamma=0.01, cost=10, probability=TRUE)) system.time(svmpol <- svm(Xfinal, Y, kernel="polynomial", degree=3, coef0=0.1, gamma=0.01, cost=10, probability=TRUE)) system.time(svmpol <- svm(Xfinal, Y, kernel="polynomial", degree=4, coef0=0.1, gamma=0.01, cost=10, probability=TRUE)) system.time(svmpol <- svm(Xfinal, Y, kernel="polynomial", degree=5, coef0=0.2, gamma=0.01, cost=10, probability=TRUE)) summary(svmpol) svmpol$rho pred.svm <- predict(svmpol, newdata=testreduced, decision.values=T) # taux d'erreur en test sum(pred.svm != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.svm) # recherche des meilleurs paramètres du SVM svmpol <- tune.svm(Xfinal, Y, kernel="polynomial", degree=3, gamma=10^(-3:0), cost=10^(-1:2), coef0=seq(0,1,by=0.1), validation.x = testreduced, validation.y = factor(test[,1]), tunecontrol = tune.control(sampling = "fix", fix=1)) svmpol$best.model summary(svmpol) system.time(svmpol <- tune.svm(Xfinal, Y, kernel="polynomial", degree=3, gamma=10^(-3:-1), cost=10^(0:2), coef0=seq(0,0.2,by=0.1), validation.x = testreduced, validation.y = factor(test[,1]), tunecontrol = tune.control(sampling = "fix", fix=1))) svmpol$best.model summary(svmpol) # détermination des coefficients des prédicteurs - avec LiblineaR - noyau linéaire install.packages('LiblineaR') library(LiblineaR) s <- scale(Xfinal) dim(Y) m <- LiblineaR(data=s, target=Y, type=4, cost=0.1, epsilon = 0.01, bias=TRUE, verbose=TRUE) # application du score SVM #st <- scale(testreduced, apply(Xfinal,2,mean), apply(Xfinal,2,sd)) st <- scale(testreduced, attr(s,"scaled:center"), attr(s,"scaled:scale")) pred.svm <- predict(m, newx=st, decisionValues=TRUE) # taux d'erreur en test sum(pred.svm$predictions != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.svm$predictions) # --------------------------------------------------------------------------------------------------------- # Chapitre 17 : Reconnaissance de chiffres manuscrits avec un réseau de neurones (base MNIST) # --------------------------------------------------------------------------------------------------------- train <- read.csv("C:/Users/.../mnist_train.csv", header=TRUE) test <- read.csv("C:/Users/.../mnist_test.csv", header=TRUE) # affichage de 100 chiffres par(mfrow = c(10,10), mai = c(0,0,0,0)) for(i in 1:100) { y = as.matrix(train[i, 2:785]) dim(y) = c(28, 28) image(y[,nrow(y):1], axes = FALSE, col = gray(255:0/255)) text(0.2, 0, train[i,1], cex = 3, col = "grey", pos = c(3,4)) } # Perceptron : modélisation directe (sans passer par une ACP) # avec une variable à expliquer facteur : pas besoin de préciser une sortie softmax (nb_poids <- (784+10)*1 + (10+1)) # calcul du nb de pondérations du réseau (nb_poids <- (784+10)*2 + (10+2)) # calcul du nb de pondérations du réseau library(nnet) set.seed(235) #rn <- nnet(label ~ ., data=train, size=1, maxit=100) system.time(rn <- nnet(factor(label) ~ ., data=train, size=1, maxit=100)) system.time(rn <- nnet(factor(label) ~ ., data=train, size=1, decay=1, maxit=100)) system.time(rn <- nnet(factor(label) ~ ., data=train, size=2, maxit=100, MaxNWts =10000)) system.time(rn <- nnet(factor(label) ~ ., data=train, size=4, maxit=100, MaxNWts =10000, decay=0.5)) system.time(rn <- nnet(factor(label) ~ ., data=train, size=5, maxit=100, MaxNWts =10000, decay=1)) # d'après Hastie-Tibshirani-Friedman, "The Elements of Statistical Learning", 2d ed. page 407, on atteint 20% de taux d'erreur avec une couche cachée de 2570 unités pred.rn <- predict(rn, newdata=test[,-1], type="class") table(pred.rn) # taux d'erreur en test sum(pred.rn != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.rn) # Perceptron sur facteurs ACP avec sortie softmax set.seed(235) trainDataset <- train[sample(nrow(train)),] # mélange des images X <- trainDataset[,-1] Y <- trainDataset[,1] Xreduced <- X/255 # ACP Xcov <- cov(Xreduced) pcaX <- prcomp(Xcov) str(pcaX) summary(pcaX) dim(pcaX$rotation) # valeurs propres plot(pcaX) screeplot(pcaX) # les deux plots sont identiques ev <- pcaX$sdev^2 # valeurs propres head(ev) head(eigen(cov(Xcov))$values) head(ev/sum(ev)) head(cumsum(ev/sum(ev)), 10) # sommes cumulées des variances plot(ev[1:100]) # matrices d'apprentissage en entrée du réseau de neurones nbaxes <- 58 dim(as.matrix(Xreduced)) # matrice avec les 60000 images en lignes et les 784 pixels en colonnes dim(pcaX$rotation[,1:nbaxes]) # matrice avec les 784 pixels en lignes les 58 axes en colonnes Xfinal <- as.matrix(Xreduced) %*% pcaX$rotation[,1:nbaxes] # les 44 premières composantes représentent 99% de la variance # les 58 premières composantes représentent 99,5% de la variance # les 82 premières composantes représentent 99,8% de la variance dim(Xfinal) class(Y) length(Y) head(Y) Y <- class.ind(Y) # le vecteur des 60000 chiffres est remplacé par une matrice de 60000 lignes et 10 colonnes indicatrices, une par chiffre, # avec une valeur 0 dans chaque colonne, sauf pour la colonne correspondant au chiffre qui aura la valeur 1 rm(Xreduced) gc() # données de test transformées par ACP testreduced <- test[,-1]/255 testreduced <- as.matrix(testreduced) %*% pcaX$rotation[,1:nbaxes] # Perceptron sur axes factoriels (nb_poids <- (nbaxes+10)*23 + (10+23)) # calcul du nb de pondérations du réseau (nb_poids <- (nbaxes+10)*400 + (10+400)) # calcul du nb de pondérations du réseau set.seed(235) #system.time(rn <- nnet(Xfinal, Y, size=23, softmax=TRUE, maxit=100, MaxNWts =10000)) #system.time(rn <- nnet(Xfinal, Y, size=100, softmax=TRUE, maxit=100, MaxNWts =10000, decay=0)) system.time(rn <- nnet(Xfinal, Y, size=400, softmax=TRUE, maxit=300, MaxNWts =30000, decay=1)) system.time(rn <- nnet(Xfinal, Y, size=400, softmax=TRUE, maxit=800, MaxNWts =80000)) system.time(rn <- nnet(Xfinal, Y, size=600, softmax=TRUE, maxit=1000, MaxNWts =50000)) system.time(rn <- nnet(Xfinal, Y, size=23, softmax=TRUE, maxit=100, MaxNWts =10000)) system.time(rn <- nnet(Xfinal, Y, size=23, softmax=TRUE, maxit=100, MaxNWts =10000, decay=1)) # decay = 0, 1, 0.5, 0.1, 0.01, 0.001, 0.0001 rn pred.rn <- predict(rn, newdata=testreduced, type="class") # taux d'erreur en test sum(pred.rn != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.rn) > system.time(rn <- nnet(Xfinal, Y, size=100, softmax=TRUE, maxit=300, MaxNWts =10000, decay=1)) # weights: 6910 initial value 330578.531491 iter 10 value 56654.927274 iter 20 value 22084.007444 iter 30 value 13419.245075 iter 40 value 10058.478566 iter 50 value 8465.091930 iter 60 value 7592.456396 iter 70 value 7012.319682 iter 80 value 6498.125377 iter 90 value 6020.497475 iter 100 value 5601.653886 iter 110 value 5269.171002 iter 120 value 5044.433003 iter 130 value 4883.719878 iter 140 value 4759.164540 iter 150 value 4667.044692 iter 160 value 4610.378902 iter 170 value 4550.123460 iter 180 value 4505.435596 iter 190 value 4426.000959 iter 200 value 4373.691262 iter 210 value 4303.246784 iter 220 value 4211.707418 iter 230 value 4148.896972 iter 240 value 4080.020313 iter 250 value 4029.580392 iter 260 value 3978.909042 iter 270 value 3932.588902 iter 280 value 3896.150698 iter 290 value 3863.673253 iter 300 value 3828.594977 final value 3828.594977 stopped after 300 iterations utilisateur système écoulé 2476.34 0.40 2490.40 > rn a 58-100-10 network with 6910 weights options were - softmax modelling decay=1 > pred.rn <- predict(rn, newdata=testreduced, type="class") > # taux d'erreur en test > sum(pred.rn != test[,1]) / nrow(test) [1] 0.0182 > # matrice de confusion > table(test[,1], pred.rn) pred.rn 0 1 2 3 4 5 6 7 8 9 0 972 0 1 0 1 2 1 1 1 1 1 0 1127 1 1 0 1 2 2 1 0 2 6 0 1009 3 1 0 5 6 2 0 3 0 0 1 997 0 2 0 4 2 4 4 1 1 3 0 964 0 1 0 1 11 5 2 1 0 10 0 875 3 0 0 1 6 5 3 2 1 4 2 938 0 3 0 7 0 5 10 2 2 0 0 1005 1 3 8 4 0 2 5 2 1 3 2 952 3 9 2 4 0 3 9 4 1 5 2 979 > system.time(rn <- nnet(Xfinal, Y, size=200, softmax=TRUE, maxit=300, MaxNWts =15000, decay=1)) # weights: 13810 initial value 352099.042087 iter 10 value 45300.882111 iter 20 value 23213.645471 iter 30 value 13380.257098 iter 40 value 9334.920319 iter 50 value 7298.442633 iter 60 value 6176.065736 iter 70 value 5539.131982 iter 80 value 5069.304874 iter 90 value 4699.493606 iter 100 value 4379.703268 iter 110 value 4119.569757 iter 120 value 3923.645568 iter 130 value 3781.416833 iter 140 value 3671.144142 iter 150 value 3585.016457 iter 160 value 3517.275723 iter 170 value 3466.152179 iter 180 value 3421.854444 iter 190 value 3360.567472 iter 200 value 3308.918806 iter 210 value 3276.838832 iter 220 value 3234.333307 iter 230 value 3206.473782 iter 240 value 3189.380050 iter 250 value 3166.000449 iter 260 value 3133.798630 iter 270 value 3102.622382 iter 280 value 3080.720897 iter 290 value 3065.716845 iter 300 value 3054.716661 final value 3054.716661 stopped after 300 iterations utilisateur système écoulé 5800.46 1.05 5830.61 > rn a 58-200-10 network with 13810 weights options were - softmax modelling decay=1 > pred.rn <- predict(rn, newdata=testreduced, type="class") > # taux d'erreur en test > sum(pred.rn != test[,1]) / nrow(test) [1] 0.0155 > # matrice de confusion > table(test[,1], pred.rn) pred.rn 0 1 2 3 4 5 6 7 8 9 0 973 0 1 0 0 0 2 1 3 0 1 0 1129 1 0 0 1 2 0 2 0 2 3 0 1012 1 3 0 2 7 4 0 3 0 0 3 997 0 2 0 4 1 3 4 1 1 1 0 966 0 4 1 1 7 5 3 0 0 6 0 876 3 1 2 1 6 3 2 0 1 6 4 940 0 2 0 7 0 3 7 1 0 0 0 1012 2 3 8 4 1 2 4 1 1 1 1 956 3 9 3 2 0 3 7 4 1 5 0 984 > > system.time(rn <- nnet(Xfinal, Y, size=300, softmax=TRUE, maxit=300, MaxNWts =30000, decay=1)) # weights: 20710 initial value 406255.722777 iter 10 value 58754.980989 iter 20 value 19359.529540 iter 30 value 12038.566529 iter 40 value 8617.528527 iter 50 value 6895.734833 iter 60 value 5954.079694 iter 70 value 5243.446666 iter 80 value 4660.912646 iter 90 value 4246.016184 iter 100 value 3941.970757 iter 110 value 3721.805478 iter 120 value 3561.879128 iter 130 value 3450.610284 iter 140 value 3370.009658 iter 150 value 3308.614455 iter 160 value 3262.168329 iter 170 value 3215.018591 iter 180 value 3183.704335 iter 190 value 3150.744116 iter 200 value 3122.228875 iter 210 value 3097.351746 iter 220 value 3079.625871 iter 230 value 3065.927469 iter 240 value 3052.641062 iter 250 value 3035.093812 iter 260 value 3017.282154 iter 270 value 2999.052553 iter 280 value 2986.148288 iter 290 value 2977.322731 iter 300 value 2969.623986 final value 2969.623986 stopped after 300 iterations utilisateur système écoulé 11832.99 2.26 11935.49 > rn a 58-300-10 network with 20710 weights options were - softmax modelling decay=1 > pred.rn <- predict(rn, newdata=testreduced, type="class") > # taux d'erreur en test > sum(pred.rn != test[,1]) / nrow(test) [1] 0.014 > # matrice de confusion > table(test[,1], pred.rn) pred.rn 0 1 2 3 4 5 6 7 8 9 0 975 0 1 0 0 1 1 1 1 0 1 0 1127 2 1 0 0 2 1 2 0 2 3 0 1011 2 2 0 2 6 6 0 3 0 0 3 997 0 3 0 2 4 1 4 1 1 2 0 969 0 2 1 0 6 5 2 0 0 8 0 875 3 1 3 0 6 3 2 0 1 3 2 946 0 1 0 7 0 3 6 1 0 0 0 1013 1 4 8 2 1 2 2 0 1 1 1 962 2 9 2 2 0 4 7 6 0 3 0 985 > > system.time(rn <- nnet(Xfinal, Y, size=400, softmax=TRUE, maxit=300, MaxNWts =30000, decay=1)) # weights: 27610 initial value 564383.877939 iter 10 value 67677.721852 iter 20 value 22857.651146 iter 30 value 13641.065716 iter 40 value 9443.754462 iter 50 value 7666.397670 iter 60 value 6739.849288 iter 70 value 5956.022808 iter 80 value 5137.962240 iter 90 value 4488.416830 iter 100 value 4077.325874 iter 110 value 3802.631008 iter 120 value 3628.680390 iter 130 value 3517.441152 iter 140 value 3443.991116 iter 150 value 3386.781148 iter 160 value 3341.926154 iter 170 value 3310.880243 iter 180 value 3280.382326 iter 190 value 3252.384293 iter 200 value 3225.390599 iter 210 value 3199.626354 iter 220 value 3174.102526 iter 230 value 3156.589028 iter 240 value 3138.488050 iter 250 value 3126.367897 iter 260 value 3113.461026 iter 270 value 3097.190611 iter 280 value 3077.375781 iter 290 value 3062.464888 iter 300 value 3051.628529 final value 3051.628529 stopped after 300 iterations utilisateur système écoulé 18900.51 5.89 19010.41 > rn a 58-400-10 network with 27610 weights options were - softmax modelling decay=1 > pred.rn <- predict(rn, newdata=testreduced, type="class") > # taux d'erreur en test > sum(pred.rn != test[,1]) / nrow(test) [1] 0.0144 > # matrice de confusion > table(test[,1], pred.rn) pred.rn 0 1 2 3 4 5 6 7 8 9 0 974 0 1 0 0 0 2 1 1 1 1 0 1130 2 1 0 0 1 0 1 0 2 3 0 1013 2 1 0 2 6 5 0 3 0 0 2 998 0 4 0 2 3 1 4 1 0 1 0 966 0 4 0 1 9 5 2 0 0 6 0 877 3 1 3 0 6 3 2 0 1 3 4 943 0 2 0 7 0 4 7 1 0 0 0 1009 2 5 8 2 0 2 3 1 1 2 2 958 3 9 2 2 0 2 7 4 1 3 0 988 > # affichage des 100 premières erreurs de classement err100 <- head(which(pred.rn != test[,1]), 100) par(mfrow = c(10,10), mai = c(0,0,0,0)) for(i in err100) { y = as.matrix(test[i, 2:785]) dim(y) = c(28, 28) image(y[,nrow(y):1], axes = FALSE, col = gray(255:0/255)) text(0.2, 0, test[i,1], cex = 3, col = gray(0.75), font=3, pos = c(3,4)) text(0.2, 1, pred.rn[i], cex = 3, col = gray(0.25), pos = c(1,2)) } # --------------------------------------------------------------------------------------------------------- # Chapitre 18 : Un exemple de credit scoring avec PYTHON # --------------------------------------------------------------------------------------------------------- # Import des librairies import pandas as pd import numpy as np import matplotlib.pyplot as plt import os import time # Lecture des données url = "https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data" columns = [ "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" ] df = pd.read_csv(url, sep=" ", names=columns) df.Cible = df.Cible - 1 # pour avoir une cible entre 0 et 1 df.drop(columns=["Etranger"], inplace=True) # Dictionnaires de recodage pour chaque colonne recodes = { 'Comptes': { 'A11': 'CC débiteur', 'A12': 'Solde entre 0 et 200 euros', 'A13': 'Solde plus de 200 euros', 'A14': 'Pas de compte' }, 'Historique_credit': { 'A30': 'Impayés passés', 'A31': 'Impayé en cours dans autre banque', 'A32': 'Jamais aucun crédit', 'A33': 'Crédits en cours sans retard', 'A34': 'Crédits passés sans retard' }, 'Objet_credit': { 'A40': 'Voiture neuve', 'A41': 'Voiture occasion', 'A42': 'Mobilier', 'A43': 'Vidéo HIFI', 'A44': 'Electroménager', 'A45': 'Travaux', 'A46': 'Etudes', 'A47': 'Vacances', 'A48': 'Formation', 'A49': 'Business', 'A410': 'Autres' }, 'Epargne': { 'A61': 'Moins de 100 euros', 'A62': 'Entre 100 et 500 euros', 'A63': 'Entre 500 et 1000 euros', 'A64': 'Plus de 1000 euros', 'A65': 'Sans épargne' }, 'Anciennete_emploi': { 'A71': 'Sans emploi', 'A72': 'Empl moins de 1 an', 'A73': 'Empl entre 1 et 4 ans', 'A74': 'Empl dans 4 et 7 ans', 'A75': 'Empl plus de 7 ans' }, 'Situation_familiale': { 'A91': 'Homme divorcé-séparé', 'A92': 'Femme divorcée-séparée-mariée', 'A93': 'Homme célibataire', 'A94': 'Homme marié-veuf', 'A95': 'Femme célibataire' }, 'Garanties': { 'A101': 'Pas de garantie', 'A102': 'Co-emprunteur', 'A103': 'Garant' }, 'Biens': { 'A121': 'Immobilier', 'A122': 'Assurance-vie', 'A123': 'Voiture ou autre', 'A124': 'Aucun bien connu' }, 'Autres_credits': { 'A141': 'Autres banques', 'A142': 'Établissements crédit', 'A143': 'Aucun crédit' }, 'Statut_domicile': { 'A151': 'Locataire', 'A152': 'Propriétaire', 'A153': 'Logement gratuit' }, 'Type_emploi': { 'A171': 'Sans emploi', 'A172': 'Non qualifié', 'A173': 'Employé-Ouvrier qualifié', 'A174': 'Cadre' }, 'Anciennete_domicile': { 1: 'Moins de 1 an', 2: 'Entre 1 et 4 ans', 3: 'Entre 4 et 7 ans', 4: 'Plus de 7 ans' }, 'Nb_credits': { 1: '1 crédit', 2: '2 ou 3 crédits', 3: '4 ou 5 crédits', 4: 'Plus de 6 crédits' }, 'Nb_pers_charge': { 1: '0 ou 2 pers', 2: 'Plus de 3 pers' }, 'Taux_effort': { 1: 'Moins de 20 pct', 2: 'Entre 20 et 25 pct', 3: 'Entre 25 et 35 pct', 4: 'Plus de 35 pct' }, 'Telephone': { 'A191': 'Sans Tél', 'A192': 'Avec Tél' } } # application des recodages au DataFrame for col, mapping in recodes.items(): df[col] = df[col].replace(mapping) # affichage du résultat print(df) # Nombre de valeurs manquantes par variable missing_values = df.isna().sum() print(f"Nombre de valeurs manquantes :\n{missing_values}\n") # Histogrammes des variables explicatives # création d'une figure matplotlib avec autant de sous-graphiques que de colonnes dans le dataframe import math nbcol = 4 fig, axs = plt.subplots(int(math.ceil(len(df.columns) / nbcol)), nbcol, figsize=(20, 30)) for i, col in enumerate(df.columns): axs[i // nbcol, i % nbcol].hist(df[col], bins=10, color='lightgray') axs[i // nbcol, i % nbcol].set_title(col) axs[i // nbcol, i % nbcol].tick_params(axis='x', rotation=90, labelsize=12) # rotation des étiquettes de l'axe des abscisses plt.tight_layout(rect=[0, 0.03, 1, 0.95]) plt.show() # transformation des variables qualitatives => format "category" (y compris la variable à expliquer) for col in df.columns: if col not in ['Duree_credit', 'Montant_credit', 'Age']: #df[col] = pd.Categorical(df[col], categories=df[col].unique()) df[col] = df[col].astype('category') # affichage des taux de défaut par modalités - valeurs d'information import pandas as pd import matplotlib.pyplot as plt def calculate_iv(df, feature, target): """Calcule l'Information Value (IV) pour une variable catégorielle.""" df_woe = df.groupby(feature, observed=True)[target].agg(['sum', 'count']) df_woe['non_event'] = df_woe['count'] - df_woe['sum'] df_woe['f_non_event'] = df_woe['non_event'] / df_woe['non_event'].sum() df_woe['f_event'] = df_woe['sum'] / df_woe['sum'].sum() df_woe['woe'] = np.log(df_woe['f_event'] / df_woe['f_non_event']) df_woe['iv'] = (df_woe['f_event'] - df_woe['f_non_event']) * df_woe['woe'] return df_woe['iv'].sum() # sélection des colonnes catégorielles autre que "Cible" categorical_columns = [col for col in df.select_dtypes(include='category').columns if col != 'Cible'] # définition de la disposition des sous-graphiques : 4 colonnes et autant de lignes que nécessaire n_cols = 4 n_rows = (len(categorical_columns) + n_cols - 1) // n_cols # Calcul du nombre de lignes fig, axes = plt.subplots(n_rows, n_cols, figsize=(15, n_rows * 5), constrained_layout=True) # aplatissement de la grille d'axes pour l'itération axes = axes.flatten() for idx, col in enumerate(categorical_columns): # calcul de la proportion de `Cible = 1` pour chaque modalité proportions = df[df['Cible'] == 1].groupby(col, observed=True).size() / df.groupby(col, observed=False).size() # calcul de l'IV iv = 100*calculate_iv(df.assign(Cible=df['Cible'].astype(int)), col, 'Cible') # Convert 'Cible' to int within calculate_iv # histogramme sur l'axe spécifié ax = axes[idx] proportions.plot(kind='bar', color='gray', edgecolor='black', ax=ax) ax.axhline((df['Cible'] == 1).mean(), color='red', linestyle='--') # Ligne horizontale en pointillé à 0.3 ax.set_xlabel(col) ax.set_ylabel("Proportion de Cible = 1") ax.set_ylim(0, 0.7) # Limite des proportions entre 0 et 1 ax.tick_params(axis='x', rotation=90) # ajout de la valeur d'information dans le coin supérieur droit ax.text(0.95, 0.95, f"VI = {iv:.2f}", transform=ax.transAxes, ha='right', va='top', bbox=dict(facecolor='white', edgecolor='none', alpha=0.7)) # suppression des axes inutilisés si le nombre de graphiques est inférieur à la grille for i in range(len(categorical_columns), len(axes)): fig.delaxes(axes[i]) plt.show() # Liaisons entre les variables explicatives quantitatives : diagramme de paires import seaborn as sns sns.pairplot(df[['Duree_credit', 'Montant_credit', 'Age']].join(df['Cible']), hue='Cible') # Liaisons entre les variables explicatives qualitatives : V de Cramer import matplotlib.pyplot as plt import numpy as np from scipy.stats import chi2_contingency # calcul des V de Cramer entre les variables catégorielles df_quali = df.select_dtypes(exclude=['number']).drop(columns=['Cible']) n = len(df_quali) V = np.zeros((len(df_quali.columns), len(df_quali.columns))) for i in range(len(df_quali.columns)): for j in range(i+1, len(df_quali.columns)): #if df_quali.iloc[:, i].dtype != "object" or df_quali.iloc[:, j].dtype != "object": # continue contingency_table = pd.crosstab(df_quali[df_quali.columns[i]], df_quali[df_quali.columns[j]]) V[i][j] = np.sqrt(chi2_contingency(contingency_table)[0] / n) V[j][i] = V[i][j] V = V.round(2) print("Taille de V:", V.shape) # création d'une matrice de corrélation (combinant les coefficients de corrélation et le V de Cramer) corr_matrix = V.copy() # affichage du graphique en couleurs représentant les corrélations entre les variables plt.figure(figsize=(16, 12)) sns.heatmap(corr_matrix, annot=True, cmap="Greys_r", center=0, fmt=".2f", xticklabels=df_quali.columns, yticklabels=df_quali.columns) plt.show() # Conversion des variables catégorielles en variables dummies (indicatrices numériques binaires) df_dummy = pd.get_dummies(df_quali).astype(int) # Xa = vecteur des indicatrices des variables catégorielles (variables continues non discrétisées) Xa = pd.concat([df[['Duree_credit', 'Montant_credit', 'Age']], df_dummy], axis=1) # mMdélisation avec LightGBM import lightgbm as lgb print('Version de lightgbm :', lgb.__version__) # Import des librairies pour le scoring from sklearn.model_selection import train_test_split from sklearn.metrics import roc_auc_score, roc_curve, auc, confusion_matrix, recall_score, precision_score, f1_score, accuracy_score # Conversion des DataFrames pandas directement en objets pour LightGBM y_train = y_train.astype('int') y_test = y_test.astype('int') # remarque : contrairement à XGBoost, lightGBM nécessite des variables à expliquer numériques (non catégorielles) lgb_train = lgb.Dataset(X_train, label=y_train) lgb_valid = lgb.Dataset(X_test, label=y_test, reference=lgb_train) # Paramètres LightGBM params = { 'device': 'cpu', # cpu ou gpu ou cuda 'objective': 'binary', 'metric': 'auc', 'learning_rate': 0.61, 'feature_fraction_bynode': 0.3, 'lambda_l1': 0.00021, 'lambda_l2': 0.17, 'num_leaves': 128, 'max_depth': 4, 'verbosity': -1, 'random_state': 123, 'num_threads': 2 # nombre de cœurs } # Entraînement du modèle evals={} start_time = time.time() lgb_model = lgb.train(params, lgb_train, valid_sets=[lgb_valid], num_boost_round=1000, callbacks=[lgb.early_stopping(stopping_rounds=10, verbose=True), lgb.record_evaluation(evals)]) print(f"Temps d'exécution: {time.time() - start_time:.2f} secondes") # Historique de l'évolution de l'AUC from lightgbm import plot_metric plot_metric(evals) # Prédictions # utilise la meilleure itération et pas la dernière y_pred = lgb_model.predict(X_test) # Calcul de la courbe ROC from sklearn.metrics import roc_curve, auc fpr, tpr, thresholds = roc_curve(y_test, y_pred) # Calcul de l'aire sous la courbe ROC (AUC) auc_value = auc(fpr, tpr) # Affichage de la courbe ROC plt.figure(figsize=(8, 6)) plt.plot(fpr, tpr, color='darkorange', lw=2, label='Courbe ROC (AUC = %0.3f)' % auc_value) plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--') plt.xlim([0.0, 1.0]) plt.ylim([0.0, 1.05]) plt.xlabel('Taux de Faux Positifs (FPR)') plt.ylabel('Taux de Vrais Positifs (TPR)') plt.title('Courbe ROC') plt.legend(loc="lower right") plt.show() # Importance des variables : liste lgm_imp = lgb_model.feature_importance(importance_type='gain') feature_names = lgb_model.feature_name() importance_df = pd.DataFrame({'Feature': feature_names, 'Gain': lgm_imp}) importance_df = importance_df.sort_values('Gain', ascending=False).head(30) # Ajout d'un numéro de ligne en première colonne importance_df = importance_df.reset_index(drop=True) importance_df.index = importance_df.index + 1 print(importance_df) # Importance des variables : graphique (utilisant matplotlib) import matplotlib.pyplot as plt plt.figure(figsize=(12, 6)) plt.bar(importance_df['Feature'], importance_df['Gain']) plt.title('Importance des variables') plt.xlabel('Variables') plt.ylabel('Gain') plt.xticks(rotation=90) plt.tight_layout() plt.show() # SHAP values for model interpretation : avec Explainer import shap explainer = shap.TreeExplainer(lgb_model) # Calcul des valeurs SHAP shap_values = explainer(X_train) # Graphique des valeurs SHAP shap.summary_plot(shap_values, features=X_train, feature_names=X_train.columns, max_display=12) # Affichage des valeurs SHAP pour un individu row_idx = 0 shap.plots.bar(shap_values[row_idx]) # --------------------------------------------------------------------------------------------------------- # Chapitre 20 : Text mining # --------------------------------------------------------------------------------------------------------- # documentation du package quanteda : # https://cran.r-project.org/web/packages/quanteda/vignettes/quickstart.html # le code qui suit a été testé avec la version 4.3.1 du package quanteda # certaines fonctions mentionnées dans le livre (metadoc, metacorpus, texts) ont disparu de quanteda à partir de sa version 3 # plus de détails ici : https://blog.quanteda.org/2021/04/07/whats-new-in-quanteda-version-3.0/ install.packages('quanteda') library(quanteda) install.packages('readtext') library(readtext) # source pour Esope : téléchargement sur https://fr.wikisource.org/wiki/Fables_d%E2%80%99%C3%89sope esope <- readtext("C:/Users/.../esope.txt", encoding="UTF-8") fontaine <- readtext("C:/Users/.../lafontaine.txt", encoding="UTF-8") # corpus composé d'un seul texte contenant toutes les fables moncorpus1 <- corpus(esope) moncorpus2 <- corpus(fontaine) summary(moncorpus1) summary(moncorpus2) moncorpus <- moncorpus1 + moncorpus2 summary(moncorpus, showmeta = TRUE) # découpage du corpus en textes : un texte par fable (elles sont séparées chacune par 4 retours à la ligne) install.packages('tokenizers') library(tokenizers) tok <- tokenize_paragraphs(moncorpus, paragraph_break = "\n\n\n\n", simplify = TRUE) # nouveau corpus composé d'un texte par fable #corp <- corpus(tok) # si un seul texte dans le corpus initial corp <- corpus(unlist(tok)) # si plusieurs textes dans le corpus initial summary(corp) corp[294] # extraction d'un texte du corpus # extraction des titres des fables (séparés par 2 espaces du reste du texte) titres <- sapply(1:ndoc(corp), function(x) substr(tokenize_paragraphs(corp[x], paragraph_break = " ")[[1]][1], 1, 60)) titres # ajout de métadonnées docvars(corp, "language") <- "french" #docvars(corp, "fable") <- paste("Fable", 1:ndoc(corp), sep = "_") docvars(corp, "chapitre") <- titres[1:ndoc(corp)] summary(corp, showmeta = TRUE) # affichage des 100 premiers documents du corpus summary(corp, 294, showmeta = TRUE) # tokenisation tokens(corp) tokens(corp, remove_punct = TRUE) # bigrammes (sans signes de ponctuation toks <- tokens_ngrams(tokens(corp, remove_punct = TRUE), n = 2, concatenator = " ") # équivalent : toks <- tokens(corp, remove_punct = TRUE) |> tokens_ngrams(n = 2, concatenator = " ") # passage en bas de casse toks <- tokens_tolower(tokens(corp)) # collocation : trouve les n-grammes qui séparent le mieux les documents (plus grand écart entre nij et mij # où nij = # occurrences observées du terme j dans le document i # et mij = # occurrences attendues du terme j dans le document i # https://quanteda.io/reference/textstat_collocations.html install.packages('quanteda.textstats') library(quanteda.textstats) textstat_collocations(toks, size = 2, min_count = 20) # mots dans leur contexte kwic(tokens(corp), "renard") # matrice documents-termes dtm <- dfm(tokens(corp)) # version synthétique dtm <- dfm(tokens_remove(tokens_tolower(tokens(corp, remove_punct = TRUE)), pattern = c("a", "à", "qu’il", "d’un", stopwords("fr")))) # version équivalente mais peut-être plus lisible dtm <- dfm( tokens(corp, remove_punct = TRUE) |> tokens_tolower() |> tokens_remove(pattern = c("a", "à", "qu’il", "d’un", stopwords("fr"))) ) dtm dtm[1:10, 1:10] # termes les plus fréquents (avec package quanteda) topfeatures(dtm,10) # associations de termes (corrélation >= 0.2) (avec package quanteda) sim <- textstat_simil(dtm, dtm[,c("loup", "renard")], method = "cosine", margin = "features") sim <- textstat_simil(dtm, dtm[,c("loup", "renard")], method = "jaccard", margin = "features") sim <- textstat_simil(dtm, dtm[,c("loup", "renard")], method = "correlation", margin = "features") lapply(as.list(sim), head, 6) # similarités entre textes distance <- textstat_simil(dtm, margin = "documents", method = "cosine") distance # diminution de la sparsité dtm presDfm <- dfm_trim(dtm, min_termfreq = 10, min_docfreq = 3) # mots apparaissant au moins 10 fois et dans au moins 3 documents presDfm presDfm[1:20, 1:10] # hierarchical clustering - get distances on normalized dfm presDistMat <- textstat_dist(dfm_weight(presDfm, "prop")) presDistMat <- textstat_dist(dfm_tfidf(presDfm)) presDistMat <- textstat_simil(dfm_tfidf(presDfm), method = "cosine") # hiarchical clustering the distance object presCluster <- hclust(presDistMat) # label with document names presCluster$labels <- docnames(presDfm) presCluster$labels <- presDfm@docvars[,2] # plot as a dendrogram plot(presCluster, xlab = "", sub = "") plot(presCluster, xlab = "", sub = "", main = "Euclidean Distance on Normalized Token Frequency") # pondération TF-IDF tfidfDtm <- dfm_tfidf(dtm) tfidfDtm[1:10,1:10] # le package text2vec contient une fonction pour créer une matrice TF-IDF # https://cran.r-project.org/web/packages/text2vec/vignettes/text-vectorization.html install.packages('text2vec') library(text2vec) dtm_tfidf <- fit_transform(dtm, TfIdf$new()) dtm_tfidf[1:10,1:10] # AFC presDfm <- dfm_trim(dtm, min_termfreq = 10, min_docfreq = 3) # mots apparaissant au moins 10 fois et dans au moins 3 documents presDfm install.packages('quanteda.textmodels') library(quanteda.textmodels) afc <- textmodel_ca(presDfm) afc <- textmodel_ca(dfm_tfidf(presDfm)) afc <- textmodel_ca(dfm_weight(presDfm, "prop")) class(afc) library(ca) plot.ca(afc) # avec FactoMineR library(FactoMineR) AFC <- CA(convert(presDfm, to = "data.frame")[, -1]) # extraction du nom des auteurs substr(row.names(dtm),1,5) strsplit(row.names(dtm), ".txt") unlist(strsplit(row.names(dtm), ".txt"))[seq(1, 2*nrow(dtm), by = 2)] docvars(corp, "auteur") <- unlist(strsplit(row.names(dtm), ".txt"))[seq(1, 2*nrow(dtm), by = 2)] summary(corp, showmeta = TRUE) # analyse factorielle multiple duale presDfm <- dfm_trim(dtm, min_termfreq = 50, min_docfreq = 10) df <- cbind(auteur=as.factor(unlist(strsplit(row.names(dtm), ".txt"))[seq(1, 2*nrow(dtm), by = 2)]), convert(presDfm, to = "data.frame")[, -1]) dim(df) head(df) # la fonction DMFA suppose que chaque variable ait une variance non nulle dans chacun des groupes défini par la variable "auteur" (num.fact = 1) # sinon : Erreur dans eigen(Cov[[j]]) : valeurs infinies ou manquantes dans 'x' (cov() génère une matrice singulière → valeurs propres infinies) AFMD <- DMFA(df, num.fact=1, ncp=2) # il faut alors supprimer toutes les colonnes qui sont constantes dans au moins un groupe vars <- colnames(df)[-1] badcols <- character(0) for (g in levels(df[[1]])) { sub <- df[df[[1]]==g, vars, drop = FALSE] if (nrow(sub) >= 2) { v <- sapply(sub, function(x) var(as.numeric(as.character(x)), na.rm=TRUE)) badcols <- union(badcols, names(v)[!is.finite(v) | v == 0]) } else { cat("Modalité", g, "a", nrow(sub), "obs -> vérifier/traiter\n") } } cat("Colonnes constantes dans au moins un groupe:", paste(badcols, collapse = ", "), "\n") df_clean <- df[, !(names(df) %in% badcols)] dim(df_clean) AFMD <- DMFA(df_clean, num.fact=1, ncp=2) # wordcloud install.packages('quanteda.textplots') library(quanteda.textplots) set.seed(235) # seulement les termes de fréquence >= 6 # random.order = FALSE => les mots fréquents sont affichés d'abord, donc au centre textplot_wordcloud(dtm, min.freq = 6, random.order = FALSE, rot.per = .25, colors = RColorBrewer::brewer.pal(8,"Dark2")) textplot_wordcloud(dtm, min.freq = 30, random.order = FALSE) # affichage d'un "comparison.cloud" comparant les fréquences de termes à travers plusieurs documents # matrice documents-termes agrégeant les documents selon un critère dtm <- dfm(tokens_remove(tokens_tolower(tokens(corp, remove_punct = TRUE)), pattern = c("a", "à", "qu’il", "d’un", stopwords("fr")))) dtm.gp <- dfm_group(dtm, groups = auteur) dtm.gp textplot_wordcloud(dtm.gp, min.freq = 30, random.order = FALSE, comparison = TRUE, max.words=100, colors=gray.colors(3)) textplot_wordcloud(dtm.gp, min_count = 30, random_order = FALSE, comparison = TRUE, max_words=100, color=gray.colors(3)) # classement des fables (classifieur bayésien naïf) library(quanteda.textmodels) dtm <- dfm(tokens_remove(tokens_tolower(tokens(corp, remove_punct = TRUE)), pattern = c("a", "à", "qu’il", "d’un", stopwords("fr")))) dtm <- dfm_trim(dtm, min_termfreq = 10, min_docfreq = 3) dtm y <- as.factor(unlist(strsplit(row.names(dtm), ".txt"))[seq(1, 2*nrow(dtm), by = 2)]) table(y) set.seed(235) s <- sample(1:294,100) (classif <- textmodel_nb(dtm[-s,], y[-s], prior = "docfreq")) (classif <- textmodel_nb(dtm[-s,], y[-s], prior = "termfreq")) (classif <- textmodel_nb(dtm[-s,], y[-s], prior = "uniform")) # matrice de confusion table(y[-s], predict(classif)) # taux de bonne classification sum(y[-s] == predict(classif)) / length(y[-s]) # matrice de confusion table(y[s], predict(classif, dtm[s,])) # taux de bonne classification sum(y[s] == predict(classif, dtm[s,])) / length(y[-s]) # autres classifieurs library(caret) train <- as.matrix(dtm[-s,]) train <- as.data.frame(train) train <- cbind(y=y[-s], train) table(train$y) valid <- as.matrix(dtm[s,]) valid <- as.data.frame(valid) valid <- cbind(y=y[s], valid) table(valid$y) # modèle fit <- train(y ~ ., data = train, method = 'bayesglm') fit <- train(y ~ ., data = train, method = 'glmboost') fit <- train(y ~ ., data = train, method = 'rpart') fit <- train(y ~ ., data = train, method = 'rf') fit <- train(y ~ ., data = train, method = 'naive_bayes') fit <- train(y ~ ., data = train, method = 'svmRadial') fit <- train(y ~ ., data = train, method = 'svmLinear') fit <- train(y ~ ., data = train, method = 'mlp') fit <- train(y ~ ., data = train, method = 'lda') # KO # matrice de confusion table(train$y, predict(fit, newdata = train)) # taux de bonne classification sum(train$y == predict(fit, newdata = train)) / length(train$y) # matrice de confusion table(valid$y, predict(fit, newdata = valid)) # taux de bonne classification sum(valid$y == predict(fit, newdata = valid)) / length(valid$y) # réseau de neurones library(nnet) summary(train) system.time(rn <- nnet(y ~ ., data=train, size=1, decay=1, maxit=100)) pred.rn <- predict(rn, newdata=valid, type="class") # taux d'erreur en test sum(pred.rn != valid[,"y"]) / nrow(valid) # matrice de confusion table(valid[,"y"], pred.rn) # classement avec le package RTextTools install.packages('RTextTools') library(RTextTools) set.seed(235) s <- sample(1:nrow(dtm)) data <- convert(dtm[s,], to = "data.frame") matrix <- create_matrix(data) container <- create_container(matrix, y[s], trainSize=1:200, testSize = 201:nrow(data), virgin=TRUE) # apprentissage SVM <- train_model(container,"SVM") GLMNET <- train_model(container,"GLMNET") SLDA <- train_model(container,"SLDA") BOOSTING <- train_model(container,"BOOSTING") BAGGING <- train_model(container,"BAGGING") RF <- train_model(container,"RF") TREE <- train_model(container,"TREE") # prédiction SVM_CLASSIFY <- classify_model(container, SVM) GLMNET_CLASSIFY <- classify_model(container, GLMNET) SLDA_CLASSIFY <- classify_model(container, SLDA) BOOSTING_CLASSIFY <- classify_model(container, BOOSTING) BAGGING_CLASSIFY <- classify_model(container, BAGGING) RF_CLASSIFY <- classify_model(container, RF) TREE_CLASSIFY <- classify_model(container, TREE) # résultats create_precisionRecallSummary(container, cbind(SVM_CLASSIFY, GLMNET_CLASSIFY, SLDA_CLASSIFY, BOOSTING_CLASSIFY, BAGGING_CLASSIFY, RF_CLASSIFY, TREE_CLASSIFY)) # ok analytics <- create_analytics(container, cbind(SVM_CLASSIFY, GLMNET_CLASSIFY, SLDA_CLASSIFY, BOOSTING_CLASSIFY, BAGGING_CLASSIFY, RF_CLASSIFY, TREE_CLASSIFY)) # ko summary(analytics) # --------------------------------------------------------------------------------------------------------- # Chapitre 20 : Lemmatisation # --------------------------------------------------------------------------------------------------------- # https://edutechwiki.unige.ch/fr/Tutoriel_koRpus install.packages('koRpus') library(koRpus) # utilisation de la fonction treetag faisant appel à TreeTagger à installer d'abord # https://www.cis.lmu.de/~schmid/tools/TreeTagger/#Windows setwd("C:/Users/.../Text Mining") Encoding(corp[25]) <- "UTF-8" # langue française install.koRpus.lang("fr") library("koRpus.lang.fr") tagged.text <- treetag(char_tolower(corp[25]), format="obj", treetagger="manual", lang="fr", encoding="UTF-8", TT.options=list(path="C:/TreeTagger/", preset="fr", no.unknown=FALSE)) tagged.text <- treetag("le chat a mangé la souris", format="obj", treetagger="manual", lang="fr", TT.options=list(path="C:/TreeTagger/", preset="fr", no.unknown=FALSE)) tagged.text <- treetag("le chat a mangé la souris", format="obj", treetagger="manual", lang="fr", TT.options=list(path="C:/TreeTagger/", preset="fr", no.unknown=TRUE)) tagged.text <- treetag("le chat mange les souris et les trouve bonnes.", format="obj", treetagger="manual", lang="fr", add.desc = TRUE, TT.options=list(path="C:/TreeTagger/", preset="fr")) # Un objet de texte tagué contient 5 "slots" (propriétés en langage orienté objet): lang, desc, tokens, features et feat_list. # tokens est un data frame qui contient 8 colonnes (token, tag, lemma, lttr, etc.) et N lignes (une pour chaque token). # Dans R, pour accéder aux slots d'une classe de type S4, on peut utiliser la fonction slot(slot,obj) ou plus simplement utiliser la notation @ str(tagged.text) tagged.text summary(tagged.text) plot(tagged.text) taggedText(tagged.text)$token # mots taggedText(tagged.text)$lemma # mots lemmatisés taggedText(tagged.text)$lemma[which(taggedText(tagged.text)$tag=="NAM"|substr(taggedText(tagged.text)$tag,1,3)=="VER")] # verbes et auxiliaires # texte lemmatisé paste(taggedText(tagged.text)$lemma, collapse=" ") # texte lemmatisé (avec seulement les noms et les verbes) paste(taggedText(tagged.text)$lemma[which(taggedText(tagged.text)$tag=="NOM"|substr(taggedText(tagged.text)$tag,1,3)=="VER")], collapse=" ") # lisibilité ARI(tagged.text) # diversité MTLD(tagged.text) lex.div(tagged.text) R.ld(tagged.text) # suppression du package de l'environnement remove.packages('koRpus') # --------------------------------------------------------------------------------------------------------- # Chapitre 20 : Détection de thématiques # --------------------------------------------------------------------------------------------------------- # Allocation de Dirichlet latente (LDA) install.packages('topicmodels') library(topicmodels) lda <- LDA(convert(dtm, to = "topicmodels"), k = 4, control = list(seed = 235)) lda <- LDA(dtm, k = 4, control = list(seed = 235)) lda <- CTM(dtm, k = 4, control = list(seed = 235)) lda <- CTM(dtm, k = 4, control = list(seed = 235, var = list(tol = 10^-4), em = list(tol = 10^-3))) terms(lda,10) table(topics(lda, 1)) # calcul de la log-vraisemblance de chaque modèle : à minimiser best.model <- lapply(seq(2, 6, by = 1), function(d){LDA(dtm, d, control = list(seed = 235))}) (best.model.logLik <- as.data.frame(as.matrix(lapply(best.model,logLik)))) which.min(best.model.logLik$V1) # comparaison de plusieurs méthodes LDA k <- 4 # nombre de sujets SEED <- 235 CSC_TM <- list(VEM = LDA(dtm, k = k, control = list(seed = SEED)), VEM_fixed = LDA(dtm, k = k, control = list(estimate.alpha = FALSE, seed = SEED)), Gibbs = LDA(dtm, k = k, method = "Gibbs", control = list(seed = SEED, burnin = 1000, thin = 100, iter = 1000)), CTM = CTM(dtm, k = k, control = list(seed = SEED, var = list(tol = 10^-4), em = list(tol = 10^-3)))) # le modèle le plus différenciant est celui dont l'entropie est la plus basse : ici le modèle variationnel VEM sapply(CSC_TM, function(x) mean(apply(posterior(x)$topics, 1, function(z) - sum(z * log(z))))) # c'est lié à la plus faible valeur du paramètre de concentration du modèle variationnel sapply(CSC_TM[1:3], slot, "alpha") Topic <- topics(CSC_TM[["VEM"]], 1) # no de la thématique de chaque document table(Topic) (Terms <- terms(CSC_TM[["VEM"]], 10)) (Terms <- terms(CSC_TM[["VEM_fixed"]], 10)) (Terms <- terms(CSC_TM[["Gibbs"]], 10)) (Terms <- terms(CSC_TM[["CTM"]], 10)) # Analyse sémantique latente (LSA) library(quanteda.textmodels) lsa <- textmodel_lsa(dtm, nd=4) # valeurs singulières. head(lsa$sk) # matrice des documents dim(lsa$docs) head(lsa$docs) # matrice des termes dim(lsa$features) head(lsa$features) # premiers termes de chaque thématique. toptopics <- function(x){ top <- order(lsa$features[,x],decreasing=T) row.names(lsa$features)[head(top,20)] } sapply(1:4,toptopics) # fréquence de la première thématique de chaque document. table(max.col(textmodel_lsa(dtm, nd=4)$docs)) # --------------------------------------------------------------------------------------------------------- # Chapitre 20 : Représentation Word2Vec # --------------------------------------------------------------------------------------------------------- # https://code.google.com/archive/p/word2vec/ # http://scoms.hypotheses.org/657 # Avec H2O (sur corpus des fables d'Esope et de La Fontaine) # http://docs.h2o.ai/h2o/latest-stable/h2o-docs/data-science/word2vec.html # https://github.com/h2oai/h2o-3/blob/master/h2o-r/demos/rdemo.word2vec.craigslistjobtitles.R library(h2o) h2o.init() # Création d'un H2OFrame contenant les fables (la matrice dtm a été construite plus haut) auteurs <- as.factor(unlist(strsplit(row.names(dtm), ".txt"))[seq(1, 2*nrow(dtm), by = 2)]) textes <- as.h2o(data.frame(auteur = auteurs, fable = as.character(corp))) head(textes) nrow(textes) h2o.table(textes$auteur) # Tokenisation du texte en mots words <- h2o.tokenize(textes$fable, "\\\\W+") nrow(words) head(as.data.frame(words), 20) # Construction du modèle Word2Vec w2v.model <- h2o.word2vec(words, epochs = 10) # Synonymes d'un mot h2o.findSynonyms(w2v.model, "renard", count = 10) # Calcul des vecteurs Word2Vec plot.vecs <- h2o.transform(w2v.model, words, aggregate = "average") # Création des échantillons d'apprentissage et test data <- h2o.cbind(textes["auteur"], plot.vecs) data.split <- h2o.splitFrame(data, ratios = 0.8) # Modèle de gradient boosting gbm.model <- h2o.gbm(x = names(plot.vecs), y = "auteur", training_frame = data.split[[1]], validation_frame = data.split[[2]]) summary(gbm.model) # Prédiction predict <- function(texte, w2v, gbm) { words <- tokenize(as.character(as.h2o(texte))) texte.vec <- h2o.transform(w2v, words, aggregate_method = "average") h2o.predict(gbm, texte.vec) } print(predict("les fables qu'on attribue à Ésope", w2v.model, gbm.model)) # Arrêt de l'instance H2O h2o.shutdown() # Avec package wordVectors # https://github.com/bmschmidt/wordVectors # https://github.com/bmschmidt/wordVectors/blob/master/vignettes/introduction.Rmd # https://bookworm.benschmidt.org/posts/2015-10-25-Word-Embeddings.html library(devtools) install_github("bmschmidt/wordVectors") library(wordVectors) library(magrittr) vecteur <- "C:/Users/.../enwik9.bin" vecteur <- "C:/Users/.../frWac_non_lem_no_postag_no_phrase_500_skip_cut200.bin" #if (!file.exists("cookbook_vectors.bin")) {model = train_word2vec("cookbooks.txt","cookbook_vectors.bin",vectors=200,threads=4,window=12,iter=5,negative_samples=0)} else model = read.vectors("cookbook_vectors.bin") # lecture d'un modèle existant model <- read.vectors(vecteur) # création d'un modèle à partir d'un corpus model <- train_word2vec("cookbooks.txt", "cookbook_vectors.bin", vectors=200, threads=2, window=12, iter=5, negative_samples=0) model %>% closest_to("king") (ana <- closest_to(model, model[["king"]])) model[[c("king","queen")]] (ana <- closest_to(model, model[[c("king","queen")]])) model %>% closest_to(model[[c("king","man","woman")]], 20) closest_to(model, ~ "king" - "man" + "woman", n = 20) closest_to(model, ~ "roi" - "homme" + "femme", n = 20) closest_to(model, ~ "paris" - "france" + "italy", n = 20) model %>% closest_to(model[[c("wolf","fox","rat")]], 20) (ana <- closest_to(model, model[[c("wolf","fox","rat","donkey")]], 20)) groupe <- model[[ana$word, average=FALSE]] plot(groupe, method="pca") # clustering set.seed(235) centers = 50 clustering = kmeans(model, centers=centers, iter.max = 20) sapply(sample(1:centers,10),function(n) { names(clustering$cluster[clustering$cluster==n][1:10]) }) # Avec H2O (sur titres d'offres d'emploi) library(h2o) h2o.init() job.titles.path = "https://raw.githubusercontent.com/h2oai/sparkling-water/rel-1.6/examples/smalldata/craigslistJobTitles.csv" job.titles <- h2o.importFile(job.titles.path, destination_frame = "jobtitles", col.names = c("category", "jobtitle"), col.types = c("Enum", "String"), header = TRUE) head(as.data.frame(job.titles), 10) nrow(job.titles) h2o.table(job.titles$category) STOP_WORDS = c("ax","i","you","edu","s","t","m","subject","can","lines","re","what", "there","all","we","one","the","a","an","of","or","in","for","by","on", "but","is","in","a","not","with","as","was","if","they","are","this","and","it","have", "from","at","my","be","by","not","that","to","com","org","like","likes","so") tokenize <- function(sentences, stop.words = STOP_WORDS) { tokenized <- h2o.tokenize(sentences, "\\\\W+") # convert to lower case tokenized.lower <- h2o.tolower(tokenized) # remove short words (less than 2 characters) tokenized.lengths <- h2o.nchar(tokenized.lower) tokenized.filtered <- tokenized.lower[is.na(tokenized.lengths) || tokenized.lengths >= 2,] # remove words that contain numbers tokenized.words <- tokenized.filtered[h2o.grep("[0-9]", tokenized.filtered, invert = TRUE, output.logical = TRUE),] # remove stop words tokenized.words[is.na(tokenized.words) || (! tokenized.words %in% STOP_WORDS),] } predict <- function(job.title, w2v, gbm) { words <- tokenize(as.character(as.h2o(job.title))) job.title.vec <- h2o.transform(w2v, words, aggregate_method = "AVERAGE") h2o.predict(gbm, job.title.vec) } # tokenisation words <- tokenize(job.titles$jobtitle) words nrow(words) head(as.data.frame(words), 10) # construction du modèle word2vec w2v.model <- h2o.word2vec(words, sent_sample_rate = 0, epochs = 10) # vérification : synonymes du mot 'teacher' print(h2o.findSynonyms(w2v.model, "teacher", count = 5)) # calcul d'un vecteur pour chaque texte (= titre d'annonce) # the input is treated as sequences of words delimited by NA. # Each word of a sequences is internally mapped to a vector, and vectors belonging to the same sentence are averaged and returned in the result job.title.vecs <- h2o.transform(w2v.model, words, aggregate_method = "AVERAGE") job.title.vecs dim(job.title.vecs) # échantillonage (en supprimant les textes inconnus donc sans plongement vectoriel) valid.job.titles <- ! is.na(job.title.vecs$C1) nrow(valid.job.titles) data <- h2o.cbind(job.titles[valid.job.titles, "category"], job.title.vecs[valid.job.titles, ]) data.split <- h2o.splitFrame(data, ratios = 0.8) # modèle de gradient boosting gbm.model <- h2o.gbm(x = names(job.title.vecs), y = "category", training_frame = data.split[[1]], validation_frame = data.split[[2]]) summary(gbm.model) # exemples de prédictions print(predict("school teacher having holidays every month", w2v.model, gbm.model)) print(predict("developer with 3+ Java experience, jumping", w2v.model, gbm.model)) print(predict("Financial accountant CPA preferred", w2v.model, gbm.model)) # arrêt de l'instance H2O h2o.shutdown()