diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..232a77286ed9bf6ec9bb89b7a500cba5dc6ae7af Binary files /dev/null and b/.DS_Store differ diff --git a/Data_Song.RDataTmp b/Data_Song.RDataTmp new file mode 100755 index 0000000000000000000000000000000000000000..b51a0a2ba950dc0850bb9cefecb97dae7293e389 Binary files /dev/null and b/Data_Song.RDataTmp differ diff --git a/Description_COR_KHA.Rmd b/Description_COR_KHA.Rmd new file mode 100755 index 0000000000000000000000000000000000000000..48ea37aa49619b31dee50de3d27c3ec757abaed6 --- /dev/null +++ b/Description_COR_KHA.Rmd @@ -0,0 +1,177 @@ +--- +title: "Projet - Mathématiques Régression Régularisée - Rapport n°2" +author: "Thibault Cordier and Julien Khamphousone" +date: "12/10/2017" +output: pdf_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + + +```{r echo = F} +setwd("~/Documents/Projet/Projets-ENSIIE-MA/MRR") +rm(list=ls()) +graphics.off() + +require(ggplot2) +require(GGally) + +#lm.ridge +require(MASS) + +#glmnet +require(glmnet) +require(corrplot) +``` + +## Importation de nos données d'entrainement et de test + +```{r} +my_colnames <- rep(NA, 91) +my_colnames[1] <- "Year" +for(i in (1:12)) my_colnames[i+1] <- paste("T_Ave_",i,sep="") +for(i in (1:78)) my_colnames[i+13] <- paste("T_Cov_",i,sep="") + + +data_train <- read.table("YearPredictionMSD.txt", sep=",", nrows = 463715, col.names = my_colnames) +data_test <- read.table("YearPredictionMSD.txt", sep=",", col.names = my_colnames, skip = 463715) +``` + +## Séparation de nos données en classes + +```{r} + +f1 <- function(x,q,repr){ + if(which.max(x <= q) - 1 == 0) repr[1][1] + else repr[which.max(x <= q) - 1][1] +} +q <- quantile(data_train$Year, probs = 0:10/10) +repr <- (q[2:length(q)] + q[1:length(q)-1])/2 + +class_year <- apply(data_train["Year"], 1, f1, q, repr) + +data_train["Year"] <- class_year + +class_year <- apply(data_test["Year"], 1, f1, q, repr) + +data_test["Year"] <- class_year +``` + + +## Sélection partielle de nos données d'entrainement et de test +```{r} +n_train <- dim(data_train)[1] +n_test <- dim(data_test)[1] +sample_train <- data_train[sample(n_train,1000),] +sample_test <- data_test[sample(n_test,1000),] +``` + +## Données manquantes +```{r} +sum(is.na(data_train[,1])) +``` + + +## Visualisation des données +```{r} +ggpairs(sample_train[,c(1,sample(2:91,5))], + lower = list(continuous = wrap("points", alpha = 0.5, size=0.15))) + +ggplot(sample_train, aes(Year, T_Ave_1)) + geom_point() + ggtitle("Scatterplot : +Year en fonction de T_Ave_1") + +{ i <- floor(runif(1, 2, 91)) ; +ggplot(sample_train, aes_string("Year", my_colnames[i])) + geom_point() + + ggtitle(paste("Scatterplot : Year = f(", my_colnames[i], ")", sep="")) +ggally_box_no_facet(sample_train, aes_string("Year", my_colnames[i])) +} + +summary(data_train["Year"]) + +ggplot(data_train, aes(Year)) + geom_bar() + + ggtitle("Barplot : Distribution des chansons selon leur année de sortie") + +ggplot(data_train, aes(1,Year)) + geom_boxplot() + + ggtitle("Boxplot : Distribution des chansons selon leur année de sortie") + +reg <- lm(Year~., data_train) +plot(reg$model$Year, reg$residuals, xlab = "Year", ylab = "Résidu") +``` + +## Modèle de régression régularisé + +```{r, echo = F} +reg <- lm(Year~., sample_train) +``` + +```{r} +summary(reg) +``` + +## Modèle de régression Lasso-Ridge-Elastic + +```{r} +x.train <- as.matrix(scale(sample_train[2:length(sample_train)], center = TRUE)) +x.test <- as.matrix(scale(sample_test[2:length(sample_test)], center = T)) + +y.train <- sample_train[,c("Year")] +y.test <- sample_test[,c("Year")] +fit.lasso <- glmnet(x.train, y.train, family = "multinomial",alpha=1) +fit.ridge <- glmnet(x.train, y.train, family = "multinomial",alpha=0) +fit.elnet <- glmnet(x.train, y.train, family = "multinomial",alpha=0.5) + +for (i in 0:2) { + assign(paste("fit", i, sep=""), cv.glmnet(x.train, y.train, type.measure="mse", + alpha=i/2,family="multinomial")) +} + +par(mfrow=c(3,2)) + +plot(fit.lasso, xvar="lambda", main="lasso") +plot(fit2, main="LASSO") + +plot(fit.ridge, xvar="lambda", main="ridge") +plot(fit0, main="Ridge") + +plot(fit.elnet, xvar="lambda", main="elnet") +plot(fit1, main="Elastic Net") + +``` + + +```{r} +yhat0 <- predict(fit0, s=fit0$lambda.1se, newx=x.test) +yhat1 <- predict(fit1, s=fit1$lambda.1se, newx=x.test) +yhat2 <- predict(fit2, s=fit2$lambda.1se, newx=x.test) + +mse0 <- mean((yhat0)^2) +mse1 <- mean((yhat1)^2) +mse2 <- mean((yhat2)^2) + +mse0 ; mse1 ; mse2 +``` + +## Modèle de régression Ridge + +```{r, echo = F} +Regression_Ridge <- lm.ridge(Year~., data.frame(scale(sample_train, center = FALSE)), lambda = seq(1,100,0.1)) +``` + +```{r} +coef(Regression_Ridge) +ridge_plot = data.frame(Regression_Ridge$GCV,Regression_Ridge$lambda) +names(ridge_plot) = c("GCV","lambda") +ggplot(data = ridge_plot, aes(x = lambda, y = GCV)) + geom_point() +``` + + +```{r} +which.min(Regression_Ridge$GCV) +coefridge = coef(lm.ridge(Year~., sample_train, lambda = 687)) +coefridge +``` + + + diff --git a/Rapport 1/.DS_Store b/Rapport 1/.DS_Store new file mode 100755 index 0000000000000000000000000000000000000000..853b2c8986bb723960f5034bbe6cf319abf27341 Binary files /dev/null and b/Rapport 1/.DS_Store differ diff --git a/Rapport 1/Rapport_1.Rmd b/Rapport 1/Rapport_1.Rmd new file mode 100755 index 0000000000000000000000000000000000000000..d66a4773911634974a3765a0aa46f6b5efc59749 --- /dev/null +++ b/Rapport 1/Rapport_1.Rmd @@ -0,0 +1,109 @@ +--- +title: "Projet MRR - YearPredictionMSD Data Set" +author: "Thibault Cordier & Julien Khamphousone" +date: "29 November 2017" +output: pdf_document +geometry: "left=2cm,right=2cm,top=2cm,bottom=2cm" +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +```{r, include=FALSE} +setwd("~/Documents/Projet/Projets-ENSIIE-MA/MRR") +rm(list=ls()) +graphics.off() + +require(ggplot2) +require(GGally) + +my_colnames <- rep(NA, 91) +my_colnames[1] <- "Year" +for(i in (1:12)) my_colnames[i+1] <- paste("T_Ave_",i,sep="") +for(i in (1:78)) my_colnames[i+13] <- paste("T_Cov_",i,sep="") + +# Importation de nos données d'entrainement et de test 463715 +n_train <- dim(data_train)[1] +sample_train <- data_train[sample(n_train,1000),] +n_test <- dim(data_test)[1] +sample_test <- data_test[sample(n_test,1000),] + +# Classification +q <- quantile(data_train$Year, probs = 0:10/10) +repr <- (q[2:length(q)] + q[2:length(q)-1])/2 + +f1 <- function(x,q,repr){ repr[which.max(x <= q) - 1][1]} + +classif <- apply(data_train["Year"],1,f1, q,repr) +``` + +# Présentation de la base de données + +La base de données (BdD) sélectionnée est issu de [**Million Song Dataset**](http://labrosa.ee.columbia.edu/millionsong/) une base de données regroupant les caractéristiques de plus d'un million de chansons, issue d'une collaboration entre **LabROSA** (Columbia University) et **The Echo Nest**. + +Elle est constituée de 515.345 observations décrites selon 91 attributs réels, soit un total de 46.381.050 informations à manipuler. Ces données nous renseignent sur l'année de sortie ainsi que sur certaines caractéristiques audio relatives aux chansons. Ces chansons observées sont principalement des pistes musicales occidentales, commercialisées entre 1922 et 2011, avec une plus forte concentration de pistes issues des années 2000. + +D'après [*UCI Machine Learning Repository*](https://archive.ics.uci.edu/ml/datasets/yearpredictionmsd#), la nature de cette BdD est adaptée à l'utilisation de méthodes de régression à des fins de prédictions. L'absence de données devra aussi être pris en compte dans l'analyse des données. + +# Nature du problème + +Un son musical se distingue essentiellement du bruit par le fait qu'il est organisé. La théorie de la musique décrit les sons musicaux avec quatre caractéristiques : la durée, la hauteur, l'intensité et le timbre. + +Dans le [*cadre de la collecte des données de **The Echo Nest***](http://docs.echonest.com.s3-website-us-east-1.amazonaws.com/_static/AnalyzeDocumentation.pdf), un analyseur décompose l'audio en éléments musicalement pertinents, séquencés dans le temps appelés *segments*. Un *segment* est une entité sonore (en général d'une durée de moins d'une seconde), relativement uniforme en timbre et en harmonie. + +Ainsi, les caractéristiques audio de la BdD, qui ont été extraites des caractéristiques *timbre* de l'API **The Echo Nest**, sont une synthèse des données relatives aux *segments* d'une chanson. La BdD reprend les moyennes et les covariances de différentes caractéristiques associées au "timbre" (de haut niveau d'abstraction) sur tous les *segments* d'une chanson, chaque *segment* étant décrit par un vecteur de timbre de 12 dimensions. + +L'hypothèse de ce projet est de considérer qu'une chanson peut être reconnaisable en partie par son timbre, et que le timbre peut potentionnellement porter la marque de l'année de sortie d'une chanson. + +Le problème ici est donc mettre en place un modèle de regression permettant d'émettre des prédictions sur l'année de sortie d'une chanson à partir de la connaissance de certaines caractéristiques audio liées au timbre de la chanson. La variable cible de notre modèle de régression sera l'année de sortie d'une chanson, allant de 1922 à 2011, étant données 12 variables explicatives relatives aux moyennes et 78 relatives aux covariances du timbre de la chanson. + +# Statistiques descriptives + +Etant donné que la taille de la BdD est imposante (515.345 observations pour 91 attributs), l'analyse statistique multivariée devient difficile à réaliser. Néanmoins, l'analyse statistique sur la variable "Year" représentant l'année de sortie d'une chanson peut nous apporter des éléments de rélexion pour la suite. + +Effectivement, la distribution des chansons selon leur année de sortie semble être non uniforme entre 1922 et 2011. Il semble même qu'elle soit déséquilibré puique les observations sont concentrées dans les années 2000. Cette distribution semble même suivre une tendance exponentielle. Cette remarques est valables pour les deux BdD. + +* 100% des chansons analysées sont sorties entre 1922 et 2011, soit sur une durée de 89 ans. +* 75 % des chansons analysées sont sorties entre 1994 et 2011, soit sur une durée de 17 ans. +* 50 % des chansons analysées sont sorties entre 2002 et 2011, soit sur une durée de 9 ans. + +```{r, echo=T} +summary(data_train["Year"]) +``` + +<!-- Ceci peut s'expliquer par le fait que la production et l'enregistrement musicale se soit intensifiée avec le temps (technologie mieux adaptée) et qu'il a été surement difficile de récupérer d'anciennes archives musicales à exploiter. Mais ce problème ne doit pas être comme une simple anecdocte lors de la construction de notre modèle de régression puisque les "vieilles" chansons seront surement considérées comme des outliers par le modèle. --> + +# Objectif du projet (Problématiques ?) + +L’objectif de ce projet est de construire un modèle prédictif utilisant des techniques de régression sur des données réelles (une méthode de régression utilisée présentée en cours et une méthode non introduite en cours à choisir parmi un ensemble de méthodes proposées). + +Dans notre cas, l'objectif sera donc d'émettre des prédictions sur l'année de sortie d'une chanson connaissant ses caractéristiques audio liées au timbre. Au cours de la mise en place de notre modèle prédictif, nous allons surement rencontrer différents problèmes. Il faudra notamment penser à : + +* prendre en compte la taille de la BdD pour des raisons de performances en terme de complexité temporelle +* prendre en compte le fait que les années de sortie des chansons soient concentrées autour des années 2000, ce qui pourrait faire révéler des outliers à tord (une chanson sorti en 1922 n'est pas une abération) +* prendre en compte l'ordre de grandeur de l'erreur réalisée par notre modèle lors de la prédiction et fixer un seuil de telle sorte qu'une certaine marge d'erreur soit acceptable ou non + +# Construction d'un premier modèle + +Le premier modèle prédictif construit a été obtenu en utilisant une régression linéraire sur notre BdD. Par une analyse de l'erreur de prédiction réalisée en fonction de l'année de sortie d'une chanson, nous remarquons que le résidu est croissant en fonction de l'année. + +Cette observation vient confirmer notre crainte concernant la considération des "vieilles" chansons comme des outliers par le modèle. Il faudra donc songer à l'avenir à mettre en place un modèle qui prenne compte de ce problème, soit par l'utilisation d'une pondération astucieuse sur les observations, soit en utilisant une méthode de régression adaptée à ce problème. + +```{r, echo=T, fig.height=3} +ggplot(data_train, aes(Year)) + geom_bar() + + ggtitle("Barplot : Distribution des chansons selon leur année de sortie") +``` +```{r, echo=T, fig.height=3} +ggplot(data_train, aes(1,Year)) + geom_boxplot() + + ggtitle("Boxplot : Distribution des chansons selon leur année de sortie") +``` + +```{r, echo=T} +reg <- lm(Year~., data_train) +``` + +```{r, echo=T, fig.height=3} +plot(reg$model$Year, reg$residuals, xlab = "Year", ylab = "Résidu") +``` + diff --git a/Rapport 1/Rapport_1.pdf b/Rapport 1/Rapport_1.pdf new file mode 100755 index 0000000000000000000000000000000000000000..b646913891fbf645ce34011c98a2e31a9f5efc42 Binary files /dev/null and b/Rapport 1/Rapport_1.pdf differ diff --git a/Rapport 2/.DS_Store b/Rapport 2/.DS_Store new file mode 100755 index 0000000000000000000000000000000000000000..9a6e234f4c404ab3ad608defb400c1a1752e4dd6 Binary files /dev/null and b/Rapport 2/.DS_Store differ diff --git a/Rapport 2/Rapport_2.Rmd b/Rapport 2/Rapport_2.Rmd new file mode 100755 index 0000000000000000000000000000000000000000..9da8c75f291fe53c309a3a758a04e8d2bdbf7dd2 --- /dev/null +++ b/Rapport 2/Rapport_2.Rmd @@ -0,0 +1,408 @@ +--- +title: "Projet MRR - YearPredictionMSD Data Set - Rapport N°2" +author: "Thibault Cordier and Julien Khamphousone" +date: "21 Décembre 2017" +output: pdf_document +geometry: "left=2cm,right=2cm,top=2cm,bottom=2cm" +--- + +```{r setup, include=F} +knitr::opts_chunk$set(echo = TRUE) +``` + +<!-- ## Initialisation des données & Chargement des packages --> + +```{r, include=FALSE} +setwd("~/Documents/Projet/Projets-ENSIIE-MA/MRR") +rm(list=ls()) +graphics.off() + +require(ggplot2) +require(GGally) + +#lm.ridge +require(MASS) + +#glmnet +require(glmnet) +require(corrplot) + +# Renommages des variables de nos données +my_colnames <- rep(NA, 91) +my_colnames[1] <- "Year" +for(i in (1:12)) my_colnames[i+1] <- paste("T_Ave_",i,sep="") +for(i in (1:78)) my_colnames[i+13] <- paste("T_Cov_",i,sep="") + +# Importation de nos données d'entrainement et de test +data_train <- read.table("YearPredictionMSD.txt", sep=",", nrows = 463715, col.names = my_colnames) +data_test <- read.table("YearPredictionMSD.txt", sep=",", col.names = my_colnames, skip = 463715) + +n_train <- dim(data_train)[1] +n_test <- dim(data_test)[1] + +sample_train <- data_train[sample(n_train,40000),] +sample_test <- data_test[sample(n_test,30000),] + +``` + +# Objectif du projet + +L’objectif de ce projet est de construire un modèle prédictif utilisant des techniques de régression sur des données réelles (une méthode de régression utilisée présentée en cours et une méthode non introduite en cours à choisir parmi un ensemble de méthodes proposées). + +Dans notre cas, l’objectif est d’émettre des prédictions sur l’année de sortie d’une chanson connaissant ses caractéristiques audio liées au timbre. Pour cela, nous nous proposons de définir deux sous ensembles de données, un d'entrainement et un de test, selon les recommandations des auteurs de la base de données. + +Avant de mettre en place nos modèles prédictifs, afin de résoudre les problématiques évoquées dans la précédente itération (taille de la base de données, distribution asymétrique de nos observations, etc.), nous avons décidé d'effectuer une sélection de nos données d'entrainement. + +# Discussion sur une première sélection de nos données d'entrainement + +La visualisation de la distribution des chansons selon leur année de sortie montre qu'elle est fortement hétérogène et dissymétrique. De plus, la visualisation de cette distribution par une boîte à moustaches nous indique aussi la présence d'un grand nombre d'outliers. (**Nous conseillons aux lecteurs de se référer au rapport précédent pour mieux comprendre cette analyse**). De plus, le modèle linéaire étudié précedemment est selon nous non suffisant puisque les résidus sont croissants en fonction de l'année. + +Suite à ces remarques, et supposant que cela est dû aux outliers nous avons envisagé deux choix de sélection des observations : + +* Soit une répartition des observations sous 10 classes de périodes (10 classes représentées par une année qui est le centre de la classe) telles qu'elles aient chacune le même nombre d'observations +* Soit une sélection homogène des observations par année telle que chaque année contienne le même nombre d'observations et en ne gardant que les années comprises entre 1970 et 2010 (afin de ne plus avoir d'outliers) + +```{r, include=FALSE} +# Répartition sous 10 périodes +### Préparation à l'étiquettage "Year" +etiq_year <- data_train["Year"] + +f1 <- function(x,q,repr){ + if(which.max(x <= q) - 1 == 0) repr[1][1] + else repr[which.max(x <= q) - 1][1] +} +q <- quantile(data_train$Year, probs = 0:10/10) +repr <- (q[2:length(q)] + q[1:length(q)-1])/2 + +class_year <- apply(data_train["Year"], 1, f1, q, repr) + +data_train["Year"] <- class_year + +sample_train_10_classes <- data_train[sample(n_train,40000),] + +data_train["Year"] <- etiq_year +``` + +```{r, include=FALSE} +# Sélection homogène sur chaque année + +# Sélection au préalable des observations (outliers) +#1994 - 1.5*(2006-1994) # Méthode des moustaches : 1976 +# 1970 -> 2010 + +selec_obs <- function(db,n) { + new_db <- NULL + for(y in 1970:2010) { + indk <- which(db[,"Year"]==y) + new_db <- rbind(new_db, db[sample(indk,min(n,length(indk))),]) + } + new_db +} + +sample_train_by_year <- selec_obs(data_train,1000) +sample_test_by_year <- selec_obs(data_test,1000) +``` + +Dans ce rapport, notre choix s'est porté sur la deuxième méthode (sélection homogène par année entre 1970 et 2010) puisque cette sélection d'observations ne changent en rien nos objectifs de prédiction que nous avons défini. Nous avons donc chercher à savoir si cette méthode nous permettait dans un premier temps d'obtenir de meilleurs résultats. + +# Modèle de régression linéaire : Analyse des performances selon la sélection des observations + +Nous avons réalisé un modèle de régression linéaire sur un échantillon aléatoire ainsi que sur un échantillon aléatoire homogène entre 1970 et 2010. + +Pour comparer les performances entre les deux modèles, nous avons décidé d'utiliser quatre outils de comparaisons des performances : + +* la visualisation par les résidus en fonction de l'année de sortie +* la visualisation par comparaison entre l'observation et la prédiction de la variable expliquée (Y et hatY) +* le calcul du MSE (Mean Square Error) sur les bases de données d'entrainement et de test +* la méthode des k-fold + +Bien que les deux modèles semblent avoir visuellement les mêmes performances en termes de prédictions (ou qu'il n'apparaît pas de différences significatives), l'erreur moyenne sur la prédiction de l'année dans le modèle avec sélection des observations est un peu plus élevée que dans le modèle sans sélection ; et ceci quelque soit la nature de la base de données de test, que ce soit avec un échantillon aléatoire hétérogène ou aléatoire homogène entre 1970 et 2010. + +```{r, include=FALSE} +# Modèles +res_all <- glm(Year~., data=sample_train) +res_all_by_year <- glm(Year~., data=sample_train_by_year) + +# Prédictions +predict_all <- predict.glm(res_all, newdata = sample_test) +predict_all_by_year <- predict.glm(res_all_by_year,newdata = sample_test) +predict_all_2 <- predict.glm(res_all, newdata = sample_test_by_year) +predict_all_by_year_2 <- predict.glm(res_all_by_year,newdata = sample_test_by_year) +``` + +```{r, echo=T, fig.height=3} +# Résidus +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(res_all$model$Year, res_all$residuals, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Train", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(res_all_by_year$model$Year, res_all_by_year$residuals, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Train", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +plot(sample_test$Year, sample_test$Year - predict_all, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Test1", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(sample_test$Year, sample_test$Year - predict_all_by_year, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Test1", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +plot(sample_test_by_year$Year, sample_test_by_year$Year - predict_all_2, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Test2", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(sample_test_by_year$Year, sample_test_by_year$Year - predict_all_by_year_2, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Test2", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +# Y par rapport à Y^ +par(mfrow=c(2,2),mar = rep(2, 4)) +plot(sample_test$Year, predict_all, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Sans sélection ; Test1", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test$Year, predict_all_by_year, + xlab = "Y", ylab = "Y^ by year", main="Y/Y^ ; Avec sélection ; Test1", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test_by_year$Year, predict_all_2, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Sans sélection ; Test2", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test_by_year$Year, predict_all_by_year_2, + xlab = "Y", ylab = "Y^ by year", main="Y/Y^ ; Avec sélection ; Test2", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +``` + +```{r, include=FALSE} +# MSE +mse_all <- sqrt(mean((res_all$model$Year - res_all$fitted.values)^2)) +mse_year <- sqrt(mean((res_all_by_year$model$Year - res_all_by_year$fitted.values)^2)) + +mse_all_p1 <- sqrt(mean((predict_all - sample_test$Year)^2)) +mse_year_p1 <- sqrt(mean((predict_all_by_year - sample_test$Year)^2)) + +mse_all_p2 <- sqrt(mean((predict_all_2 - sample_test_by_year$Year)^2)) +mse_year_p2 <- sqrt(mean((predict_all_by_year_2 - sample_test_by_year$Year)^2)) + +# KFold +kfold <- function(tab,var_exp,k) { + n <- nrow(tab) + m <- n %/% k + X <- tab[sample(nrow(tab)),] + error_test <- c() + error_train <- c() + E_train <- 0 + E_test <- 0 + + for(i in 1:k){ + ind_sample = ((i-1)*m+1):(i*m) + mod = glm(as.formula(paste(var_exp,"~.")), data=X[-ind_sample,]) + + pY = predict.glm(mod, X[-ind_sample,]) + E_train = E_train + sqrt(mean((X[-ind_sample,var_exp] - pY)^2)) + error_train <- c(error_train,sqrt(mean((X[-ind_sample,var_exp] - pY)^2))) + + pY = predict.glm(mod, X[ind_sample,]) + E_test = E_test + sqrt(mean((X[ind_sample,var_exp] - pY)^2)) + error_test <- c(error_test,sqrt(mean((X[ind_sample,var_exp] - pY)^2))) + } + E_train <- E_train/k + E_test <- E_test/k + + res <- list("Etrain" = E_train, "Etest" = E_test, + "ErrorTrain" = error_train, "ErrorTest" = error_test, + "n" = n, "m" = m, "k" = k) + return(res) +} + +mse_kf <- kfold(sample_train, "Year", 10) +mse_kf_by_year <- kfold(sample_train_by_year, "Year", 10) +#mse_kf +#mse_kf_by_year +``` + +Nous avons donc privilégié le choix de ne pas réaliser de sélection sur les observations en amont puisque une sélection de cette nature ne permet pas d'avoir de meilleures prédictions de manière significative. + +Pour poursuivre, nous avons chercher à analyser les performances de différents modèles après sélection de variables et ceci sans sélection d'observations. + +# Modèle de régression linéaire avec sélection de variables de type "forward" + +Nous avons réalisé un modèle de régression linéaire sur un échantillon aléatoire avec sélection de variables de type "forward" car nous suspectons certaines variables d'être corrélées entre elles, ce qui serait certainement la cause de la croissance des erreurs résiduelles selon l'année. + +Pour comparer les performances entre ce modèle et celui de référence, nous avons décidé d'utiliser les trois premiers outils de comparaisons des performances que précédemment (la méthode k-fold étant trop gourmande en complexité temporelle). + +**A titre indicatif, la construction du modèle de régression linéaire avec la méthode "forward" a pris environ 30 miniutes sur nos ordinateurs personnels.** + +Et étonnement, les deux modèles présentent les mêmes performances en termes de prédictions. Les graphes de visualisations des performances ainsi que l'erreur moyenne sur la prédiction sont presque surement identiques. De cette remarque nous pouvons conclure que les variables de notre base de données ne sont pas fortement corrélées et que ceci n'explique pas le problème des résidus. + +```{r, include=FALSE} +# res_all <- glm(Year~., data=sample_train) +# res_nul <- glm(Year~1, data=sample_train) +# fw_all <- step(res_nul, list(upper=res_all), direction="forward") +# +# predict_all <- predict.glm(res_all, newdata = sample_test) +# predict_all_fw <- predict.glm(fw_all, newdata = sample_test) +# +# # Résidus +# par(mfrow=c(2,2),mar = rep(2, 4)) +# plot(res_all$model$Year, res_all$residuals, +# xlab = "Y", ylab = "résidus", main="Résidus ; Standard ; Train", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# plot(fw_all$model$Year, fw_all$residuals, +# xlab = "Y", ylab = "résidus by year", main="Résidus ; Forward ; Train", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# +# plot(sample_test$Year, sample_test$Year - predict_all, +# xlab = "Y", ylab = "résidus", main="Résidus ; Standard ; Test", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# plot(sample_test$Year, sample_test$Year - predict_all_fw, +# xlab = "Y", ylab = "résidus", main="Résidus ; Forward ; Test", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# +# # Y par rapport à Y^ +# par(mfrow=c(1,2),mar = rep(2, 4)) +# plot(sample_test$Year, predict_all, +# xlab = "Y", ylab = "Y^", main="Y/Y^ ; Standard ; Test", cex = 0.3, +# xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +# plot(sample_test$Year, predict_all_fw, +# xlab = "Y", ylab = "Y^", main="Y/Y^ ; Forward ; Test", cex = 0.3, +# xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +# +# # MSE +# sqrt(mean((res_all$model$Year - res_all$fitted.values)^2)) # 9.624977 +# sqrt(mean((fw_all$model$Year - fw_all$fitted.values)^2)) # 9.626477 +# +# sqrt(mean((predict_all - sample_test$Year)^2)) # 9.45553 +# sqrt(mean((predict_all_fw - sample_test$Year)^2)) # 9.457589 +# +# formula(fw_all) +# Year ~ T_Ave_1 + T_Cov_2 + T_Ave_3 + T_Ave_2 + T_Cov_26 + T_Ave_6 + +# T_Cov_8 + T_Cov_11 + T_Cov_45 + T_Cov_15 + T_Cov_1 + T_Cov_24 + +# T_Cov_29 + T_Cov_28 + T_Ave_9 + T_Cov_47 + T_Ave_11 + T_Cov_36 + +# T_Cov_4 + T_Cov_7 + T_Cov_66 + T_Cov_73 + T_Cov_13 + T_Cov_52 + +# T_Cov_3 + T_Ave_8 + T_Cov_6 + T_Cov_46 + T_Ave_5 + T_Cov_76 + +# T_Cov_63 + T_Cov_59 + T_Cov_41 + T_Cov_53 + T_Cov_64 + T_Cov_57 + +# T_Cov_12 + T_Cov_21 + T_Cov_27 + T_Cov_72 + T_Cov_77 + T_Cov_38 + +# T_Cov_62 + T_Cov_35 + T_Cov_30 + T_Cov_19 + T_Cov_33 + T_Cov_75 + +# T_Cov_20 + T_Cov_56 + T_Ave_10 + T_Cov_70 + T_Cov_51 + T_Cov_34 + +# T_Cov_58 + T_Cov_61 + T_Cov_71 + T_Cov_50 + T_Cov_48 + T_Cov_60 + +# T_Cov_18 + T_Ave_4 + T_Cov_14 + T_Cov_40 + T_Cov_49 +``` + +Nous avons donc privilégié le choix de ne pas réaliser "systématiquement" de sélection de type "forward" sur les variables puisque une sélection de cette nature permet d'obtenir des performances semblables au cas standard. + +Mais il suffira de réduire le nombre de variables ; et de ne garder que celles significatives ; dans nos modèles pour garder les mêmes performances tout en diminuant le nombre de calculs. + +Pour poursuivre, nous avons chercher à analyser les performances de différents modèles après sélection de variables de type Lasso, Ridge et ElasticNet, avec les mêmes objectifs que précédemment. + +# Modèle de régression linéaire avec sélection de variables de type "Lasso-Ridge-ElasticNet" + +Les méthodes de type "Lasso-Ridge-ElasticNet" sont des méthodes de contraction des coefficients dans un modèle de régression linéaire sous contraintes. La méthode "Elastic-Net" a été introduite pour surmonter les limitations du "Lasso" en combinant les contraintes "Lasso" et "Ridge". + +Nous avons donc réalisé 5 modèles de régression linéaire avec sélection de variables de type "Lasso-Ridge-ElasticNet" en donnant plus ou moins d'importance sur les contraintes "Lasso" et "Ridge" ; dont 2 modèles ayant pour seule contrainte "Lasso" ou "Ridge". + +Pour comparer les performances entre ces modèles, nous avons décidé d'utiliser les mêmes outils que précédemment. + +```{r, include=FALSE} +x.train <- as.matrix(scale(sample_train[2:length(sample_train)])) +x.test <- as.matrix(scale(sample_test[2:length(sample_test)])) + +y.train <- sample_train[,c("Year")] +y.test <- sample_test[,c("Year")] + +for (i in 0:4) { + assign(paste("fit", i, sep=""), cv.glmnet(x.train, y.train, type.measure="mse", + alpha=i/4, family="gaussian")) + print(i); +} + +#fit0$lambda.min ; fit1$lambda.min ; fit2$lambda.min ; fit3$lambda.min ; fit4$lambda.min + +#plot(fit4, main="LASSO") +#plot(fit3, main="Elastic Net 0.75") +#plot(fit2, main="Elastic Net 0.5") +#plot(fit1, main="Elastic Net 0.25") +#plot(fit0, main="Ridge") + +# Modèles +fit.ridge <- glmnet(x.train, y.train, family = "gaussian", alpha=0, lambda = fit0$lambda.min) +fit.elnet0.25 <- glmnet(x.train, y.train, family = "gaussian", alpha=0.25, lambda = fit1$lambda.min) +fit.elnet0.5 <- glmnet(x.train, y.train, family = "gaussian", alpha=0.5, lambda = fit2$lambda.min) +fit.elnet0.75 <- glmnet(x.train, y.train, family = "gaussian", alpha=0.75, lambda = fit3$lambda.min) +fit.lasso <- glmnet(x.train, y.train, family = "gaussian", alpha=1, lambda = fit4$lambda.min) + +# Prédictions +yhatall <- predict.glm(res_all, newdata = sample_test) +yhat0 <- predict(fit.ridge, s=fit0$lambda.min, newx=x.test) +yhat1 <- predict(fit.elnet0.25, s=fit1$lambda.min, newx=x.test) +yhat2 <- predict(fit.elnet0.5, s=fit2$lambda.min, newx=x.test) +yhat3 <- predict(fit.elnet0.75, s=fit3$lambda.min, newx=x.test) +yhat4 <- predict(fit.lasso, s=fit4$lambda.min, newx=x.test) +``` + +```{r, echo=T, fig.height=3} +# Résidus +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(y.test, y.test - yhatall, + xlab = "Y", ylab = "résidus", main="Résidus ; Standard ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat0, + xlab = "Y", ylab = "résidus", main="Résidus ; Ridge ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat1, + xlab = "Y", ylab = "résidus", main="Résidus ; ElNet.0,25 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat2, + xlab = "Y", ylab = "résidus", main="Résidus ; ElNet.0,50 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat3, + xlab = "Y", ylab = "résidus", main="Résidus ; ElNet.0,75 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat4, + xlab = "Y", ylab = "résidus", main="Résidus ; Lasso ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +# Y par rapport à Y^ +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(y.test, yhatall, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Standard ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat0, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Ridge ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat1, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; ElNet.0,25 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat2, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; ElNet.0,50 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat3, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; ElNet.0,75 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat4, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Lasso ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +``` + +```{r, include=FALSE} +# MSE +mse_all <- sqrt(mean((y.test-yhatall)^2)) +mse_ridge <- sqrt(mean((y.test-yhat0)^2)) +mse_elastic0.25 <- sqrt(mean((y.test-yhat1)^2)) +mse_elastic0.5 <- sqrt(mean((y.test-yhat2)^2)) +mse_elastic0.75 <- sqrt(mean((y.test-yhat3)^2)) +mse_lasso <- sqrt(mean((y.test-yhat4)^2)) + +# mse_all ; mse_ridge ; mse_elastic0.25 ; mse_elastic0.5 ; mse_elastic0.75 ; mse_lasso +``` + +Et étonnement, tous les modèles présentent presque surement les mêmes performances en termes de prédictions comme précédemment. De cette remarque nous pouvons toujours conclure que les variables de notre base de données ne sont pas fortement corrélées et que ceci n'explique pas le problème des résidus. + +# Réflexions et solutions proposées pour la suite + +Nous en concluons que les modèles mis en place jusque là ne sont pas suffisants pour prédire l'année de sortie d'une chanson mieux que le modèle de régression linéaire standard. De plus, ces modèles ne répondent pas non plus au problème de résidus que nous avons rencontrés. + +Cependant, nous ne s'attendions pas à un miracle : la sélection de variables n'aurait pas pu à elle seule résoudre ce problème. Néanmoins, nous pensons de plus en plus que mettre en place un modèle de régression **linéaire** n'est pas adapté à notre base de données. + +Nous proposons pour la suite de s'orienter vers d'autres solutions pour résoudre ces problèmes persistants : +* MARS (modèles additifs et splines) +* PCA +* Modèles non linéaires + diff --git a/Rapport 2/Rapport_2.pdf b/Rapport 2/Rapport_2.pdf new file mode 100755 index 0000000000000000000000000000000000000000..355848500b8b943162bf10260f3a7dfb227b145a Binary files /dev/null and b/Rapport 2/Rapport_2.pdf differ diff --git a/Rapport Final/.DS_Store b/Rapport Final/.DS_Store new file mode 100755 index 0000000000000000000000000000000000000000..f5391679cab354b83cfc3a6e650901da1823cd8f Binary files /dev/null and b/Rapport Final/.DS_Store differ diff --git a/Rapport Final/Rapport_Final.Rmd b/Rapport Final/Rapport_Final.Rmd new file mode 100755 index 0000000000000000000000000000000000000000..0bb5250dc34c2e6e04b5faf02e9f6e3d282c3ef7 --- /dev/null +++ b/Rapport Final/Rapport_Final.Rmd @@ -0,0 +1,1288 @@ +--- +title: "Projet Mathématiques Régression Régularisée - YearPredictionMSD Data Set" +author: "Thibault Cordier & Julien Khamphousone" +date: "21 January 2018" +header-includes: + - \usepackage{tabu} + - \usepackage{multirow} +output: pdf_document + +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +```{r Initialisation Base de Données, message=FALSE, warning=FALSE, include=FALSE} +# setwd("~/Documents/Projet/Projets-ENSIIE-MA/MRR") +setwd("~/Documents/ENSIIE/Projets/Projets-ENSIIE-MA/MRR") +rm(list=ls()) +graphics.off() + +require(ggplot2) +require(GGally) +require(MASS) +require(glmnet) +require(corrplot) +require(ROCR) +require(ggfortify) + +my_colnames <- rep(NA, 91) +my_colnames[1] <- "Year" +for(i in (1:12)) my_colnames[i+1] <- paste("T_Ave_",i,sep="") +for(i in (1:78)) my_colnames[i+13] <- paste("T_Cov_",i,sep="") +data_train <- read.table("YearPredictionMSD.txt", sep=",", nrows = 463715, col.names = my_colnames) +data_test <- read.table("YearPredictionMSD.txt", sep=",", col.names = my_colnames, skip = 463715) + +n_train <- dim(data_train)[1] +n_test <- dim(data_test)[1] + +# Multiple plot function +# +# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) +# - cols: Number of columns in layout +# - layout: A matrix specifying the layout. If present, 'cols' is ignored. +# +# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), +# then plot 1 will go in the upper left, 2 will go in the upper right, and +# 3 will go all the way across the bottom. +# +multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { + require(grid) + + # Make a list from the ... arguments and plotlist + plots <- c(list(...), plotlist) + + numPlots = length(plots) + + # If layout is NULL, then use 'cols' to determine layout + if (is.null(layout)) { + # Make the panel + # ncol: Number of columns of plots + # nrow: Number of rows needed, calculated from # of cols + layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), + ncol = cols, nrow = ceiling(numPlots/cols)) + } + + if (numPlots==1) { + print(plots[[1]]) + + } else { + # Set up the page + grid.newpage() + pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) + + # Make each plot, in the correct location + for (i in 1:numPlots) { + # Get the i,j matrix positions of the regions that contain this subplot + matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) + + print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, + layout.pos.col = matchidx$col)) + } + } +} +``` + +# Présentation de la base de données + +La base de données (BdD) sélectionnée est issue de [**Million Song Dataset**](http://labrosa.ee.columbia.edu/millionsong/) une base de données regroupant les caractéristiques de plus d'un million de chansons, issue d'une collaboration entre **LabROSA** (Columbia University) et **The Echo Nest**. + +Elle est constituée de 515.345 observations décrites selon 91 attributs réels, soit un total de 46.381.050 informations à manipuler. Ces données nous renseignent sur l'année de sortie ainsi que sur certaines caractéristiques audio relatives aux chansons. Ces chansons observées sont principalement des pistes musicales occidentales, commercialisées entre 1922 et 2011, avec une plus forte concentration de pistes issues des années 2000. + +D'après [*UCI Machine Learning Repository*](https://archive.ics.uci.edu/ml/datasets/yearpredictionmsd#), la nature de cette BdD est adaptée à l'utilisation de méthodes de régressions à des fins de prédictions. + +# Nature du problème + +Un son musical se distingue essentiellement du bruit par le fait qu'il est organisé. La théorie de la musique décrit les sons musicaux avec quatre caractéristiques : la durée, la hauteur, l'intensité et le timbre. + +Dans le [*cadre de la collecte des données de **The Echo Nest***](http://docs.echonest.com.s3-website-us-east-1.amazonaws.com/_static/AnalyzeDocumentation.pdf), un analyseur décompose l'audio en éléments musicalement pertinents, séquencés dans le temps appelés *segments*. Un *segment* est une entité sonore (en général d'une durée de moins d'une seconde), relativement uniforme en timbre et en harmonie. + +Ainsi, les caractéristiques audio de la BdD, qui ont été extraites des caractéristiques *timbre* de l'API **The Echo Nest**, sont une synthèse des données relatives aux *segments* d'une chanson. La BdD reprend les moyennes et les covariances de différentes caractéristiques associées au "timbre" (de haut niveau d'abstraction) sur tous les *segments* d'une chanson, chaque *segment* étant décrit par un vecteur de timbre de 12 dimensions. + +L'hypothèse de ce projet est de considérer qu'une chanson peut être reconnaisable en partie par son timbre, et que le timbre peut potentionnellement porter la marque de l'année de sortie d'une chanson. + +Le problème ici est donc de mettre en place un modèle de regression permettant d'émettre des prédictions sur l'année de sortie d'une chanson à partir de la connaissance de certaines caractéristiques audio liées au timbre de la chanson. La variable cible de notre modèle de régression sera l'année de sortie d'une chanson, allant de 1922 à 2011, étant données 12 variables explicatives relatives aux moyennes et 78 relatives aux covariances du timbre de la chanson. + +# Statistiques descriptives + +Etant donné que la taille de la BdD est imposante (515.345 observations pour 91 attributs), l'analyse statistique multivariée devient difficile à réaliser. Néanmoins, l'analyse statistique sur la variable "Year" représentant l'année de sortie d'une chanson peut nous apporter des éléments de rélexion pour la suite. + +Effectivement, la distribution des chansons [(voir annexe 1)](#a_1) selon leur année de sortie semble être non uniforme entre 1922 et 2011. Il semble même qu'elle soit déséquilibré puique les observations sont concentrées dans les années 2000. Cette distribution semble même suivre une tendance exponentielle. Cette remarque est valables pour les deux BdD. + +* 100% des chansons analysées sont sorties entre 1922 et 2011, soit sur une durée de 89 ans. +* 75 % des chansons analysées sont sorties entre 1994 et 2011, soit sur une durée de 17 ans. +* 50 % des chansons analysées sont sorties entre 2002 et 2011, soit sur une durée de 9 ans. + +# Objectif du projet + +L’objectif de ce projet est de construire un modèle prédictif utilisant des techniques de régression sur des données réelles (une méthode de régression utilisée présentée en cours et une méthode non introduite en cours à choisir parmi un ensemble de méthodes proposées). + +Dans notre cas, l'objectif sera donc d'émettre des prédictions sur l'année de sortie d'une chanson connaissant ses caractéristiques audio liées au timbre. Au cours de la mise en place de notre modèle prédictif, nous allons surement rencontrer différents problèmes. Il faudra notamment penser à : + +* prendre en compte la taille de la BdD pour des raisons de performances en terme de complexité temporelle +* prendre en compte le fait que les années de sortie des chansons soient concentrées autour des années 2000, ce qui pourrait faire révéler des outliers à tord (une chanson sortie en 1922 n'est pas une aberration) +* prendre en compte l'ordre de grandeur de l'erreur réalisée par notre modèle lors de la prédiction et fixer un seuil de telle sorte qu'une certaine marge d'erreur soit acceptable ou non + +# Construction d'un premier modèle de régression linéaire + +Le premier modèle prédictif construit a été obtenu en utilisant une régression linéraire sur notre BdD. Le calcul de l'erreur moyenne quadratique montre que notre modèle réalise une erreur moyenne de 9.5 années, ce qui témoigne d'une performance non satisfaisante selon nous. Par une analyse de l'erreur de prédiction réalisée en fonction de l'année de sortie d'une chanson, nous remarquons aussi que le résidu est croissant en fonction de l'année [(voir annexe 2)](#a_2). + +Cette observation vient confirmer notre crainte concernant la considération des "vieilles" chansons comme des outliers par le modèle. Il faudra donc songer à l'avenir à mettre en place un modèle qui prenne compte de ce problème, soit par l'utilisation d'une pondération astucieuse sur les observations, soit en utilisant une méthode de régression adaptée à ce problème. + +Avant de mettre en place nos modèles prédictifs, afin de résoudre les problématiques évoquées précédemment, nous avons décidé d'effectuer une filtration de nos données d'entrainement. + +# Discussion sur une filtration de nos données d'entrainement + +La visualisation de la distribution des chansons selon leur année de sortie [(voir annexe 1)](#a_1) montre qu'elle est fortement hétérogène et dissymétrique. De plus, la visualisation de cette distribution par une boîte à moustaches nous indique aussi la présence d'un grand nombre d'outliers. De plus, le modèle linéaire étudié précedemment est selon nous non suffisant puisque les résidus sont croissants en fonction de l'année. + +Suite à ces remarques, et supposant que cela est dû aux outliers, nous avons envisagé deux choix de sélection des observations : + +* Soit une répartition des observations sous 10 classes de périodes (10 classes représentées par une année qui est le centre de la classe) telles qu'elles aient chacune le même nombre d'observations +* Soit une sélection homogène des observations par année telle que chaque année contienne le même nombre d'observations et en ne gardant que les années comprises entre 1970 et 2010 (afin de ne plus avoir d'outliers) + +Dans ce rapport, notre choix s'est porté sur la deuxième méthode (sélection homogène par année entre 1970 et 2010) [(voir annexe 3)](#a_3) puisque cette sélection d'observations ne change en rien nos objectifs de prédiction que nous avons défini. De plus, une sélection homogène par année entre 1970 et 2010 était raisonnable puisque nous pouvions prendre assez d'observations dans chaque année (1000 obs/year). Nous avons donc cherché à savoir si cette méthode nous permettait dans un premier temps d'obtenir de meilleurs résultats. + +# Modèle de régression linéaire : Analyse des performances selon la filtration des observations + +Nous avons réalisé un modèle de régression linéaire sur un échantillon aléatoire ainsi que sur un échantillon aléatoire homogène entre 1970 et 2010. + +Pour comparer les performances entre les deux modèles, nous avons décidé d'utiliser quatre outils de comparaisons des performances [(voir annexe 4)](#a_4) : + +* la visualisation par les résidus en fonction de l'année de sortie +* la visualisation par comparaison entre l'observation et la prédiction de la variable expliquée (Y et Y_hat) +* le calcul du MSE (Mean Square Error) sur les bases de données d'entrainement et de test +* la méthode des k-fold + +Bien que les deux modèles semblent avoir visuellement les mêmes performances en termes de prédictions (ou qu'il n'apparaît pas de différences significatives), l'erreur moyenne sur la prédiction de l'année dans le modèle avec sélection des observations est un peu plus élevée que dans le modèle sans sélection ; et ceci quelque soit la nature de la base de données de test, que ce soit avec un échantillon aléatoire hétérogène ou aléatoire homogène entre 1970 et 2010. + +Nous avons donc privilégié le choix de ne pas réaliser de sélection sur les observations en amont puisque une sélection de cette nature ne permet pas d'avoir de meilleures prédictions de manière significative. + +Pour poursuivre, nous avons chercher à analyser les performances de différents modèles après sélection de variables et sans filtration d'observations. + +# Modèle de régression linéaire avec sélection de variables "forward" + +Nous avons réalisé un modèle de régression linéaire sur un échantillon aléatoire avec sélection de variables de type *forward* car nous suspectons certaines variables d'être corrélées entre elles, ce qui serait certainement la cause de la croissance des erreurs résiduelles selon l'année. + +Pour comparer les performances entre ce modèle et celui de référence, nous avons décidé d'utiliser les trois premiers outils de comparaisons des performances que précédemment (la méthode k-fold étant trop gourmande en complexité temporelle). + +**A titre indicatif, la construction du modèle de régression linéaire avec la méthode "forward" a pris environ 30 miniutes sur nos ordinateurs personnels.** + +Et étonnement, les deux modèles présentent les mêmes performances en termes de prédictions. Les graphes de visualisations des performances ainsi que l'erreur moyenne sur la prédiction sont presque identiques. De cette remarque nous pouvons conclure que les variables de notre base de données ne sont pas fortement corrélées et que ceci n'explique pas le problème des résidus. + +Nous avons donc privilégié le choix de ne pas réaliser *systématiquement* de sélection de type *forward* sur les variables puisque une sélection de cette nature permet d'obtenir des performances semblables au cas standard. + +Mais si nous le souhaitons, il suffit de réduire le nombre de variables ; et de ne garder que celles significatives ; dans nos modèles pour garder les mêmes performances tout en diminuant le nombre de calculs. + +Pour poursuivre, nous avons chercher à analyser les performances de différents modèles après sélection de variables de type *Lasso*, *Ridge* et *ElasticNet*, avec les mêmes objectifs que précédemment. + +# Modèle de régression linéaire avec sélection de variables de type "Lasso-Ridge-ElasticNet" + +Les méthodes de type *Lasso-Ridge-ElasticNet* sont des méthodes de contraction des coefficients dans un modèle de régression linéaire sous contraintes. La méthode *Elastic-Net* a été introduite pour surmonter les limitations du *Lasso* en combinant les contraintes *Lasso* et *Ridge*. + +Nous avons donc réalisé 5 modèles de régression linéaire avec sélection de variables de type *Lasso-Ridge-ElasticNet* en donnant plus ou moins d'importance sur les contraintes *Lasso* et *Ridge* ; dont 2 modèles ayant pour seule contrainte *Lasso* ou *Ridge*. + +Pour comparer les performances entre ces modèles, nous avons décidé d'utiliser les mêmes outils que précédemment. + +Et étonnement, tous les modèles présentent presque les mêmes performances en termes de prédictions [(voir annexe 5)](#a_5) comme précédemment. De cette remarque nous pouvons toujours conclure que les variables de notre base de données ne sont pas fortement corrélées et que ceci n'explique pas le problème des résidus. + +Un résumé statistique expliquant que nos 5 modèles ont des performances similaire est l'erreur quadratique moyenne de nos modèles sur de nouvelles prédictions, résumées dans le tableau suivant : + +$$\begin{array} +{c|c|c|c|c|c} +\text{Modèle utilisé} & Ridge & E.-N.,~0.25 & E.-N.,~0.5 & E.-N.,~0.75 & Lasso\\ \hline +\text{Erreur quadratique moyenne} & 9.524821 & 9.516081 & 9.516889 & 9.512372 & 9.518247 +\end{array}$$ + +Nos modèles ont environ une erreur quadratique moyenne de 9.5, cela signifie que l'erreur moyenne sur des nouvelles prédictions est d'environ 9.5 ans, tout comme pour le modèle linéaire standard. + +# Prise en compte de la marge d'erreur sur la performance des modèles + +Nos modèles linéaires étudiés précedemment : + +* *Modèle Linéaire Généralisé* +* *Modèle Linéaire avec sélection de variables de type "forward"* +* *Modèle Ridge* +* *Modèle Lasso* +* *Modèle Elastic-Net* + +ne nous permettent pas de prédire, avec une forte performance, l'année de sortie exacte d'une musique. + +C'est pour cette raison que nous nous sommes intéressés à la prédiction, non pas sur une année de sortie exacte, mais sur une plage d'années, caractérisées par une marge d'erreur acceptable. Par exemple, pour une marge d'erreur acceptable de 4 ans, si notre musique est sortie en 2000 et que le modèle prédit 1997, nous jugerions cette prédiction satisfaisante. + +Nous récapitulons certaines marges d'erreur et le pourcentage de bonnes prédictions dans le tableau suivant, d'après le modèle de régression linéaire standard : + +$$\begin{array} +{c|c|c|c|c|c|c|c|c|c} +\text{Marge d'erreur en année} & 75 & 40 & 20 & 15 & 10 & 5 & 2 & 1 & 0,5\\ \hline +\text{Pourcentage de bonne prédiction (\%)} & 100 & 99,7 & 94,5 & 90,2 & 80,1 & 50, 0 & 21,1 & 10,6 & 5,4 +\end{array}$$ + +Nous avons aussi réalisé un graphique représantant le pourcentage de bonnes prédictions réalisées dans le modèle linéaire standard en fonction de la marge d'erreur que l'on s'accorde [(voir en annexe 6)](#a_6). Les courbes ont été obtenues sur des données d'entrainement (catégorie 1) et sur des données de test (catégore 2). + +A posteriori, accepter une marge d'erreur de 10 ans nous permet de réaliser un peu plus de 80% de bonnes prédictions sur nos données d'entrainements. Cependant, les prédictions sur de nouvelles données, forcément moins bonnes, atteignent dans notre cas un peu plus de 65%. Au vu de la complexité dont fait preuve notre BdD, ce chiffre peut nous être satisfaisant à ce stade. + +# Analyse en Composantes Principales (ACP) + +Afin de détecter quelles variables sont corrélées nous avons réalisé une Analyse en Composantes Principales. Il s'est avéré que le résultat de cette analyse est non exploitable en raison du grand nombre de variables dans notre base de données. + +En outre, projeter les individus sur le plan principal (2 variables) ou encore obtenir la projection des variables sur le cercle de corrélations n'est pas facilement interprétable car ces 2 variables ne représenteraient que 19.3% de nos données ce qui est très peu. + +En effet, le critère de Kaiser, ou critère absolu nous suggère de sélectionner toutes les valeurs propres supérieures à 1. Dans notre étude, il faudrait donc étudier 30 dimensions (ce qui représenterait alors 72.7% de l'information), ce qui n'est pas réaliste. + + +Vous trouverez néanmoins en [annexe 11](#a_11) les résultats de l'Analyse en Composantes Principales que nous avons effectuée. + +# Réflexions et solutions proposées pour la suite + +Nous en concluons que les modèles mis en place jusque là ne sont pas suffisants pour prédire l'année de sortie d'une chanson mieux que le modèle de régression linéaire standard. De plus, ces modèles ne répondent pas non plus au problème de résidus que nous avons rencontrés. + +Cependant, nous ne nous attendions pas à un miracle : la sélection de variables n'aurait pas pu à elle seule résoudre ce problème. Néanmoins, nous pensons de plus en plus que mettre en place un modèle de régression **linéaire** n'est pas adapté à notre base de données. + +Jusqu'à maintenant, l’objectif était d’émettre des prédictions sur l’année de sortie d’une chanson connaissant ses caractéristiques audio liées au timbre. Pour redéfinir nos ambitions dans ce projet, nous nous sommes demandés s'il était plutôt possible d'émettre des prédictions quant à l'appartenance d'une chanson à une période donnée, dans laquelle elle aurait pu être commercialisée. + +Ainsi, nous avons apporté de l'intérêt dans la construction de modèles prédictifs par régression logistique sur notre base de données. + +# Modèles de régression logistique : Catégoriser les chansons selon leur année de sortie + +Les modèles prédictifs construits ici ont été obtenus en utilisant une régression logistique sur notre base de données après modifications de celle-ci. Effectivement, dans le cas des régressions binomiales, la variable cible prend deux modalités possibles : {0,1}. Nous avons donc dû adapter notre base de données pour pouvoir mettre en place cette régression. + +Pour cela, il était nécessaire de fixer une année de sortie de référence *year_ref* pour catégoriser nos observations. + +* Si une chanson est sortie en une année **antérieure** à *year_ref*, alors elle recevra la modalité **0**. +* Si une chanson est sortie en une année **postérieure** à *year_ref*, alors elle recevra la modalité **1**. + +Le choix d'une année de sortie de référence s'est porté sur les trois quartiles *Q1*, *Me* et *Q3* de notre échantillon d'entrainement. En effet, ce choix nous a permis de catégoriser les observations en deux classes de modalités de taille raisonnable. + +* Pour *year_ref* = *Q1*, la modalité **0** comprend 25% des observations d'entrainement et la modalité **1** 75%. +* Pour *year_ref* = *Me*, la modalité **0** comprend 50% des observations d'entrainement et la modalité **1** 50%. +* Pour *year_ref* = *Q3*, la modalité **0** comprend 75% des observations d'entrainement et la modalité **1** 25%. + +Ainsi, chaque modalité dans chacun des cas reste *raisonnablement* représentative d'une partie des observations d'entrainement et aura donc un poids *raisonnable* dans la construction des modèles. + +Après construction des trois modèles de régression logistique, il est devenu possible de réaliser des prédictions quant à l'appartenance d'une chanson à sa période de sortie. La prédiction réalisée peut s'interpréter comme une probabilité d'appartenir à la modalité **1** contre la modalité **0**. Il est donc nécessaire de définir un seuil, arbitrairement égal à *0.5* dans un premier temps, pour décider de l'appartenance d'une chanson à une modalité. + +* Si la prédiction est **supérieure** au seuil défini, l'observation appartiendra à la modalité **1**. +* Si la prédiction est **inférieure** au seuil défini, l'observation appartiendra à la modalité **0**. + +Après réalisation des prédictions et analyse des indicateurs à travers une matrice de confusion [(voir annexe 7)](#a_7), nous constatons que les modèles possèdent des performances relativement bonnes, entre 70% et 80% de bonnes prédictions. Cependant, nous constatons aussi que le modèle à référence *Q1* présente un taux de *faux positifs* relativement à la modalité **1** élevé (supérieur à 60%) et que le modèle à référence *Q3* présente un taux de *faux négatifs* relativement à la modalité **0** élevé (supérieur à 70%). + +Nous verrons plus tard comment palier ce problème, trouver des seuils plus adaptés pour les modèles et ainsi obtenir des résultats plus satisfaisants. + +# Modèle combinant les trois modèles de regression logistique + +Les modèles prédictifs par régression logistique construits jusqu'ici nous permettent uniquement d'émettre une préference quant à l'appartenance d'une chanson parmi deux périodes données relativement à une année de référence. Par la suite, il nous a semblé intéréssant de construire un modèle regroupant ces trois modèles prédictifs. En effet, ce modèle nous permettrait de passer d'une prédiction bimodale à une prédiction multimodale (avec quatre modalités dans notre cas). + +Ainsi, avec notre nouveau modèle, nous avons donc dû adapter notre base de données de telle sorte que : + +* Si une chanson est sortie en une année **antérieure** à *Q1*, alors elle recevra la modalité **1**. +* Si une chanson est sortie en une année **postérieure** à *Q1* et **antérieure** à *Me*, alors elle recevra la modalité **2**. +* Si une chanson est sortie en une année **postérieure** à *Me* et **antérieure** à *Q3*, alors elle recevra la modalité **3**. +* Si une chanson est sortie en une année **postérieure** à *Q3*, alors elle recevra la modalité **4**. + +Pour réaliser nos prédictions sur nos observations selon notre modèle multimodal, nous avons décidé de combiner nos trois modèles logistiques déjà existants. + +* Si la prédiction réalisée par le **modèle de référence Me** est **inférieure** au seuil défini pour celui-ci, l'observation appartiendra à la modalité **1** ou **2**. + * Si la prédiction réalisée par le **modèle de référence Q1** est **inférieure** au seuil défini pour celui-ci, l'observation appartiendra à la modalité **1**. + * Sinon, l'observation appartiendra à la modalité **2**. +* Sinon, l'observation appartiendra à la modalité **3** ou **4**. + * Si la prédiction réalisée par le **modèle de référence Q3** est **inférieure** au seuil défini pour celui-ci, l'observation appartiendra à la modalité **3**. + * Sinon, l'observation appartiendra à la modalité **4**. + +Après réalisation des prédictions et analyse des indicateurs à travers une matrice de confusion [(voir annexe 8)](#a_8), nous constatons que le modèle possède une performance relativement médiocre, entre 30% et 40% de bonnes prédictions. C'est à dire que notre modèle se trompe plus d'une fois sur deux. + +Au vu de ces remarques, nous avons décidé de palier les problèmes mentionnés en cherchant des seuils plus adaptés pour nos modèles et ainsi obtenir des résultats plus satisfaisants. + +# Recherche de seuils intéressants pour la performance + +Pour nous aider dans notre choix de seuils de discrimination pertinents à la prédiction, nous avons décidé de tracer les *courbes ROC* de chacun de nos trois modèles logistiques binaires [(voir annexe 9)](#a_9). Ces courbes permettent de montrer les performances réalisées par un classificateur binaire lorsque le seuil de discrimination varie. + +Compte tenu des mauvaises performances (au sens général) de nos classificateurs, nous cherchons des seuils permettant, non pas d'obtenir une *performance* optimale dans l'absolu, mais des seuils permettant d'obtenir des taux de *faux positifs* et de *faux négatifs* relativement équilibrés et/ou relativement insensibles aux variations de la base d'observations de test et aux variations du seuil lui-même. + +Ainsi, nous avons voulu chercher le seuil associé au point de la courbe ROC ayant une dérivé égale à 1 (dans l'hypothèse d'une courbe *concave*). Effectivement : + +* Si le seuil est associé à un point de **dérivé plus grande que 1**, alors une augmentation de la valeur du seuil implique une augmentation du taux de *vrais positifs* plus élevée que celle du taux de *faux positifs*. +* Inversement, si le seuil est associé à un point de **dérivé plus petite que 1**, alors une diminution de la valeur du seuil implique une diminution du taux de *faux positifs* plus élevée que celle du taux de *vrais positifs*. + +Ainsi, tendre la valeur du seuil vers le point de la courbe ROC de dérivé 1 est pertinent dans le sens où le seuil *équilibre* le taux de *faux positifs* et de *vrais positifs* (de même que pour les *négatifs*). Concrètement, le choix de ce seuil permet de ne pas privilégier une modalité plus qu'une autre dans une certaine mesure, tout en cherchant la performance du modèle. + +Après réajustement des seuils en conséquence, réalisation des prédictions et analyse des indicateurs à travers une matrice de confusion [(voir annexe 10)](#a_10) pour les **modèles bimodaux**, nous constatons que nos nouveaux modèles possèdent des performances moins bonnes que ceux d'avant, entre 65% et 75% de bonnes prédictions. Cependant, nous constatons aussi que les taux de *faux positifs* et de *faux négatifs* sont plus équilibrés que précédemment, avec des valeurs comprises entre 20% et 35% avec moins de 10 points d'écart entre les deux taux. C'est effectivement ce que nous recherchions. + +En ce qui concerne le **modèle multimodal**, nous constatons que le nouveau modèle possède une meilleure performance que le modèle précédent, en passant de 38% à 45%. De plus, les taux de *bonnes prédictions* a priori et a posteriori sont devenus plus *équilibrés* avec ce nouveau modèle. Avec le choix de ces seuils pertinents, nous avons améliorés *globalement* les performances de notre modèle multimodal. + +# Conclusion + +L’objectif de ce projet était de construire un modèle prédictif utilisant des techniques de régression sur des données réelles. + +Il s'est avéré que la construction de modèles de régression linéaire n'était pas adaptée à notre base de données. En outre, la construction de modèles logistiques a apporté une nouvelle vision sur la base de données même si leur performances n'ont pas été à la hauteur de nos espérances. Cependant, tous ces modèles nous ont permis de comprendre la base de données et sa complexité et nous ont permis d'appréhender des méthodes permettant d'améliorer leur performances. + +\pagebreak + +# Annexes + +## Annexe 1 - Distribution des chansons selon leur année de sortie {#a_1} +```{r Distribution des chansons, echo=FALSE, fig.height=2, message=FALSE, warning=FALSE} +ggplot(data_train, aes(Year)) + geom_bar() + + ggtitle("Barplot : Distribution des chansons selon leur année de sortie") + +ggplot(data_train, aes(1,Year)) + geom_boxplot() + + ggtitle("Boxplot : Distribution des chansons selon leur année de sortie") +``` + +## Annexe 2 - Résidus pour chaque observation {#a_2} +```{r Modèle Linéaire, echo=FALSE, fig.height=2, message=FALSE, warning=FALSE} +reg <- lm(Year~., data_train) +qplot(reg$model$Year, reg$residuals, xlab = "Year", ylab = "Résidu", + main = "Scatter-plot des résidus pour chaque observation") + + geom_abline(intercept = 0, slope = 0) +``` + +## Annexe 3 - Distribution des chansons selon leur année de sortie après filtration {#a_3} +```{r Distribution après filtration, echo=FALSE, fig.height=2, message=FALSE, warning=FALSE} +n_train <- dim(data_train)[1] +n_test <- dim(data_test)[1] + +selec_obs <- function(db,n) { + new_db <- NULL + for(y in 1970:2010) { + indk <- which(db[,"Year"]==y) + new_db <- rbind(new_db, db[sample(indk,min(n,length(indk))),]) + } + new_db +} + +sample_train <- data_train[sample(n_train,50000),] +sample_test <- data_test[sample(n_test,50000),] + +sample_train_by_year <- selec_obs(data_train,1000) # sample_train contient 41000 observations +sample_test_by_year <- selec_obs(data_test,1000) + +ggplot(sample_train_by_year, aes(Year)) + geom_bar() + + ggtitle("Barplot : Distribution des chansons selon leur année de sortie après filtration") +``` + +## Annexe 4 - Modèles Prédictifs : Avec filtration {#a_4} + +```{r Modèles prédictifs avec filtration, echo=FALSE, message=FALSE, warning=FALSE} +# Modèles +res_all <- glm(Year~., data=sample_train) +res_all_by_year <- glm(Year~., data=sample_train_by_year) + +# Prédictions +predict_all <- predict.glm(res_all, newdata = sample_test) +predict_all_by_year <- predict.glm(res_all_by_year,newdata = sample_test) +predict_all_2 <- predict.glm(res_all, newdata = sample_test_by_year) +predict_all_by_year_2 <- predict.glm(res_all_by_year,newdata = sample_test_by_year) + +# Résidus +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(res_all$model$Year, res_all$residuals, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Train", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(res_all_by_year$model$Year, res_all_by_year$residuals, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Train", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +plot(sample_test$Year, sample_test$Year - predict_all, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Test1", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(sample_test$Year, sample_test$Year - predict_all_by_year, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Test1", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +plot(sample_test_by_year$Year, sample_test_by_year$Year - predict_all_2, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Test2", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(sample_test_by_year$Year, sample_test_by_year$Year - predict_all_by_year_2, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Test2", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + + +``` + +```{r Modèles prédictifs avec filtration 2, echo=FALSE, message=FALSE, warning=FALSE} +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(res_all$model$Year, res_all$fitted.values, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Sans sélection ; Train", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(res_all_by_year$model$Year, res_all_by_year$fitted.values, + xlab = "Y", ylab = "Y^ by year", main="Y/Y^ ; Avec sélection ; Train", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) + +plot(sample_test$Year, predict_all, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Sans sélection ; Test1", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test$Year, predict_all_by_year, + xlab = "Y", ylab = "Y^ by year", main="Y/Y^ ; Avec sélection ; Test1", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) + +plot(sample_test_by_year$Year, predict_all_2, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Sans sélection ; Test2", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test_by_year$Year, predict_all_by_year_2, + xlab = "Y", ylab = "Y^ by year", main="Y/Y^ ; Avec sélection ; Test2", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) + + +# MSE +mse_all <- sqrt(mean((res_all$model$Year - res_all$fitted.values)^2)) ; +cat("MSE_all : ",mse_all,"\n") +mse_year <- sqrt(mean((res_all_by_year$model$Year - res_all_by_year$fitted.values)^2)) ; +cat("MSE_year : ",mse_all,"\n") + +mse_all_p1 <- sqrt(mean((predict_all - sample_test$Year)^2)) ; +cat("MSE_all_test1 : ",mse_all_p1,"\n") +mse_year_p1 <- sqrt(mean((predict_all_by_year - sample_test$Year)^2)) ; +cat("MSE_year_test1 : ",mse_year_p1,"\n") + +mse_all_p2 <- sqrt(mean((predict_all_2 - sample_test_by_year$Year)^2)) ; +cat("MSE_all_test2 : ",mse_all_p2,"\n") +mse_year_p2 <- sqrt(mean((predict_all_by_year_2 - sample_test_by_year$Year)^2)) ; +cat("MSE_year_test2 : ",mse_year_p2,"\n") + +# KFold +kfold <- function(tab,var_exp,k) { + n <- nrow(tab) + m <- n %/% k + X <- tab[sample(nrow(tab)),] + error_test <- c() + error_train <- c() + E_train <- 0 + E_test <- 0 + + for(i in 1:k){ + ind_sample = ((i-1)*m+1):(i*m) + mod = glm(as.formula(paste(var_exp,"~.")), data=X[-ind_sample,]) + + pY = predict.glm(mod, X[-ind_sample,]) + E_train = E_train + sqrt(mean((X[-ind_sample,var_exp] - pY)^2)) + error_train <- c(error_train,sqrt(mean((X[-ind_sample,var_exp] - pY)^2))) + + pY = predict.glm(mod, X[ind_sample,]) + E_test = E_test + sqrt(mean((X[ind_sample,var_exp] - pY)^2)) + error_test <- c(error_test,sqrt(mean((X[ind_sample,var_exp] - pY)^2))) + } + E_train <- E_train/k + E_test <- E_test/k + + res <- list("Etrain" = E_train, "Etest" = E_test, + "ErrorTrain" = error_train, "ErrorTest" = error_test, + "n" = n, "m" = m, "k" = k) + return(res) +} + +mse_kf <- kfold(sample_train, "Year", 10) ; mse_kf +mse_kf_by_year <- kfold(sample_train_by_year, "Year", 10) ; mse_kf_by_year +``` + +## Annexe 5 - Modèles Prédictifs : Avec Ridge, Elastic-Net et Lasso {#a_5} +```{r Modèles prédictifs 1, echo=FALSE, message=FALSE, warning=FALSE} +sample_train <- data_train[sample(n_train,10000),] +sample_test <- data_test[sample(n_test,10000),] + +x.train <- as.matrix(scale(sample_train[2:length(sample_train)])) +x.test <- as.matrix(scale(sample_test[2:length(sample_test)])) + +y.train <- sample_train[,c("Year")] +y.test <- sample_test[,c("Year")] + +# Validation croisée +for (i in 0:4) { + assign(paste("fit", i, sep=""), cv.glmnet(x.train, y.train, type.measure="mse", + alpha=i/4,family="gaussian")) +} + +cat("Lambda min choisi par Cross-Validation, Ridge : \n", fit0$lambda.min) +cat("\nLambda min choisi par Cross-Validation, Elastic-Net 0.25 : \n", fit1$lambda.min) +cat("\nLambda min choisi par Cross-Validation, Elastic-Net 0.5 : \n", fit2$lambda.min) +cat("\nLambda min choisi par Cross-Validation, Elastic-Net 0.75 : \n", fit3$lambda.min) +cat("\nLambda min choisi par Cross-Validation, Lasso : \n", fit4$lambda.min) + +# graphics.off() +# +# autoplot(fit4, main="LASSO") +# autoplot(fit3, main="Elastic Net 0.75") +# autoplot(fit2, main="Elastic Net 0.5") +# autoplot(fit1, main="Elastic Net 0.25") +# autoplot(fit0, main="Ridge") + +res_all <- glm(Year~., data=sample_train) +fit.ridge <- glmnet(x.train, y.train, family = "gaussian",alpha=0, lambda = fit0$lambda.min) +fit.elnet0.25 <- glmnet(x.train, y.train, family = "gaussian",alpha=0.25, lambda = fit1$lambda.min) +fit.elnet0.5 <- glmnet(x.train, y.train, family = "gaussian",alpha=0.5, lambda = fit2$lambda.min) +fit.elnet0.75 <- glmnet(x.train, y.train, family = "gaussian",alpha=0.75, lambda = fit3$lambda.min) +fit.lasso <- glmnet(x.train, y.train, family = "gaussian",alpha=1, lambda = fit4$lambda.min) + +yhatall <- predict.glm(res_all, newdata = sample_test) +yhat0 <- predict(fit.ridge, s=fit0$lambda.min, newx=x.test) +yhat1 <- predict(fit.elnet0.25, s=fit1$lambda.min, newx=x.test) +yhat2 <- predict(fit.elnet0.5, s=fit2$lambda.min, newx=x.test) +yhat3 <- predict(fit.elnet0.75, s=fit3$lambda.min, newx=x.test) +yhat4 <- predict(fit.lasso, s=fit4$lambda.min, newx=x.test) + +mse_stand <- sqrt(mean((y.test-yhatall)^2)) +mse_ridge <- sqrt(mean((y.test-yhat0)^2)) +mse_elastic0.25 <- sqrt(mean((y.test-yhat1)^2)) +mse_elastic0.5 <- sqrt(mean((y.test-yhat2)^2)) +mse_elastic0.75 <- sqrt(mean((y.test-yhat3)^2)) +mse_lasso <- sqrt(mean((y.test-yhat4)^2)) + +cat("\nMean Square Error, Standard : \n", mse_stand) +cat("\nMean Square Error, Ridge : \n", mse_ridge) +cat("\nMean Square Error, Elastic-Net 0.25 : \n", mse_elastic0.25) +cat("\nMean Square Error, Elastic-Net 0.5 : \n", mse_elastic0.5) +cat("\nMean Square Error, Elastic-Net 0.75 : \n", mse_elastic0.75) +cat("\nMean Square Error, Lasso : \n", mse_lasso) +``` + +```{r Modèles prédictifs 2, echo=FALSE, message=FALSE, warning=FALSE} +# Graphiques des différents résidus +plotStandResidus <- qplot(y.test, y.test-yhatall, xlab = "Year", ylab = "Résidu", + main = "Résidus - Standard") + + geom_abline(intercept = 0, slope = 0) + +plotRidgeResidus <- qplot(y.test, y.test-yhat0, xlab = "Year", ylab = "Résidu", + main = "Résidus - Ridge") + + geom_abline(intercept = 0, slope = 0) + +plotEN025Residus <- qplot(y.test, y.test-yhat1, xlab = "Year", ylab = "Résidu", + main = "Résidus - Elastic-Net 0.25") + + geom_abline(intercept = 0, slope = 0) + +plotEN05Residus <- qplot(y.test, y.test-yhat2, xlab = "Year", ylab = "Résidu", + main = "Résidus - Elastic-Net 0.5") + + geom_abline(intercept = 0, slope = 0) + +plotEN075Residus <- qplot(y.test, y.test-yhat3, xlab = "Year", ylab = "Résidu", + main = "Résidus - Elastic-Net 0.75") + + geom_abline(intercept = 0, slope = 0) + +plotLassoResidus <- qplot(y.test, y.test-yhat4, xlab = "Year", ylab = "Résidu", + main = "Résidus - Lasso") + + geom_abline(intercept = 0, slope = 0) + +multiplot(plotStandResidus, plotRidgeResidus, plotEN025Residus, plotEN05Residus, plotEN075Residus, plotLassoResidus, cols=2) +``` + +```{r Modèles prédictifs 3, echo=FALSE, message=FALSE, warning=FALSE} + +PlotYhatStand <- qplot(y.test,yhatall, xlab = "Variable Cible", main = "Standard", + ylab = "Variable prédicte") + + geom_abline(intercept = 0, slope = 1) + +PlotYhatRidge <- qplot(y.test,yhat0, xlab = "Variable Cible", main = "Ridge", + ylab = "Variable prédicte") + + geom_abline(intercept = 0, slope = 1) + +PlotYhatEN025 <- qplot(y.test,yhat1, xlab = "Variable Cible", main = "Elastic-Net 0.25", + ylab = "Variable prédicte") + + geom_abline(intercept = 0, slope = 1) + +PlotYhatEN05 <- qplot(y.test,yhat2, xlab = "Variable Cible", main = "Elastic-Net 0.5", + ylab = "Variable prédicte") + + geom_abline(intercept = 0, slope = 1) + +PlotYhatEN075 <- qplot(y.test,yhat3, xlab = "Variable Cible", main = "Elastic-Net 0.75", + ylab = "Variable prédicte") + + geom_abline(intercept = 0, slope = 1) + +PlotYhatLasso <- qplot(y.test,yhat4, xlab = "Variable Cible", main = "Lasso", + ylab = "Variable prédicte") + + geom_abline(intercept = 0, slope = 1) + +multiplot(PlotYhatStand, PlotYhatRidge, PlotYhatEN025, PlotYhatEN05, PlotYhatEN075, PlotYhatLasso, cols=2) +``` + +## Annexe 6 - Prédictions avec marge d'erreur {#a_6} +```{r Prédictions avec marge d\'erreur, echo=FALSE, message=FALSE, warning=FALSE} +res_all <- glm(Year~., data=sample_train) +Yhat <- res_all$fitted.values +Yhat_test <- predict.glm(res_all, newdata = sample_test) +Y <- res_all$model$Year + +calcul_delta <- function(Y, Yhat, delta) { + n <- length(Y) + res <- 0 + for(i in 1:n) if(abs(Yhat[i] - Y[i]) < delta) res = res + 1 + 100*res/n +} + +MargeDerreur <- seq(0,40,0.2) +PourcentageReussiteTrain <- sapply(MargeDerreur, calcul_delta, Y=Y, Yhat=Yhat) +PourcentageReussiteTest <- sapply(MargeDerreur, calcul_delta, Y=Y, Yhat=Yhat_test) + +data_marge <- data.frame(cbind(c(MargeDerreur,MargeDerreur), + c(PourcentageReussiteTrain,PourcentageReussiteTest), + c(rep(1,length(MargeDerreur)),rep(2,length(MargeDerreur))))) +names(data_marge) <- c("Marge_Erreur","Bonnes_Predictions","Categorie") +qplot(data_marge$Marge_Erreur,data_marge$Bonnes_Predictions, color = data_marge$Categorie, + xlab = "Marge d'erreur en années", ylab = "Pourcentage de bonnes prédictions") + +``` + +## Annexe 7 - Matrice de confusion pour la logistique bimodale avec seuils arbitraires {#a_7} +```{r, fig.height=3, message=FALSE, warning=FALSE, include=FALSE} +sample_train <- data_train[sample(n_train,50000),] +sample_test <- data_test[sample(n_test,50000),] + +q_fig <- quantile(data_train$Year) + +sample_train_logic_Q1 <- sample_train +sample_train_logic_Q2 <- sample_train +sample_train_logic_Q3 <- sample_train + +sample_train_logic_Q1$Year <- ifelse(sample_train_logic_Q1$Year >= q_fig[2], 1, 0) +sample_train_logic_Q2$Year <- ifelse(sample_train_logic_Q2$Year >= q_fig[3], 1, 0) +sample_train_logic_Q3$Year <- ifelse(sample_train_logic_Q3$Year >= q_fig[4], 1, 0) + +sample_test_logic_Q1 <- sample_test +sample_test_logic_Q2 <- sample_test +sample_test_logic_Q3 <- sample_test + +sample_test_logic_Q1$Year <- ifelse(sample_test_logic_Q1$Year >= q_fig[2], 1, 0) +sample_test_logic_Q2$Year <- ifelse(sample_test_logic_Q2$Year >= q_fig[3], 1, 0) +sample_test_logic_Q3$Year <- ifelse(sample_test_logic_Q3$Year >= q_fig[4], 1, 0) + +res_logic_Q1 <- glm(Year~., data=sample_train_logic_Q1, family=binomial) +res_logic_Q2 <- glm(Year~., data=sample_train_logic_Q2, family=binomial) +res_logic_Q3 <- glm(Year~., data=sample_train_logic_Q3, family=binomial) + +# La fonction analysze.perf affiche les éléments liés à la matrice de confusion +analyze.perf <- function(res, var, seuil, data = res$model, print = 1) { + + val_ref <- t(data[var]) + val_pred <- ifelse(predict.glm(res, newdata = data, type="response")>= seuil, 1, 0) + + conf <- table(val_pred, val_ref) + + TN <- conf[1,1] ; FP <- conf[2,1] ; FN <- conf[1,2] ; TP <- conf[2,2] + av_N <- FP + TN ; av_P <- FN + TP ; ap_N <- TN + FN ; ap_P <- FP + TP + n <- TN + FP + FN + TP + + vp <- TP / av_P ; fp <- FP / av_N + vn <- TN / av_N ; fn <- FN / av_P + err <- (FN + FP) / n ; acc <- (TP + TN) / n + + if(print==1) { + cat("\n") + cat("------------------------------------------------\n") + cat("Analyse des performances du modèles\n") + cat("Comparaison entre réponses prédites et attendues\n") + cat("------------------------------------------------\n") + cat("\n Matrice de confusion :\n") ; print(conf) + cat("\n Performance : ") ; cat(acc) + cat("\n Taux d'erreur : ") ; cat(err) + cat("\n Sensibilité (vrais positifs) : ") ; cat(vp) + cat("\n Spécificité (vrais négatifs) : ") ; cat(vn) + cat("\n Taux de 'faux positifs' : ") ; cat(fp) + cat("\n Taux de 'faux négatifs' : ") ; cat(fn) + } + else{ + res <- list("conf" = conf, + "acc" = acc, + "err" = err, + "vp" = vp, + "vn" = vn, + "fp" = fp, + "fn" = fn) + return(res) + } +} + +seuil <- 0.5 +seuilQ1 <- 0.5 +seuilQ2 <- 0.5 +seuilQ3 <- 0.5 + +# Réalisation des prédictions selon un seuil arbitraire : +# Modèle bimodal sur données d'entrainement +cat("sample_train_logic : Q1 = ", q_fig[2], " | seuil = ", seuilQ1, "\n") ; analyze.perf(res_logic_Q1, "Year", seuilQ1) +cat("sample_train_logic : Me = ", q_fig[3], " | seuil = ", seuilQ2, "\n") ; analyze.perf(res_logic_Q2, "Year", seuilQ2) +cat("sample_train_logic : Q3 = ", q_fig[4], " | seuil = ", seuilQ3, "\n") ; analyze.perf(res_logic_Q3, "Year", seuilQ3) + +# Réalisation des prédictions selon un seuil arbitraire +# Modèle bimodal sur données de test +cat("sample_test_logic : Q1 = ", q_fig[2], " | seuil = ", seuilQ1, "\n") ; analyze.perf(res_logic_Q1, "Year", seuilQ1, sample_test_logic_Q1) +cat("sample_test_logic : Me = ", q_fig[3], " | seuil = ", seuilQ2, "\n") ; analyze.perf(res_logic_Q2, "Year", seuilQ2, sample_test_logic_Q2) +cat("sample_test_logic : Q3 = ", q_fig[4], " | seuil = ", seuilQ3, "\n") ; analyze.perf(res_logic_Q3, "Year", seuilQ3, sample_test_logic_Q3) + + +``` + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_train\_logic : Q1 = 1994, seuil = 0.5 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.8057 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.1943 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.9388948\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} +& \multicolumn{2}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 4146 & 2352 \\ + 1 & 7363 & 36139 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.3602398\\ +& Taux de 'faux positifs' : 0.6397602\\ +& Taux de 'faux n\'egatifs' : 0.06110519\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_train\_logic : Me = 2002, seuil = 0.5 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.72066 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.27934 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7755212\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} +& \multicolumn{2}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 15947 & 5814 \\ + 1 & 8153 & 20086 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.6617012\\ +& Taux de 'faux positifs' : 0.3382988\\ +& Taux de 'faux n\'egatifs' : 0.2244788\\ +\hline +\end{tabu}} + + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_train\_logic : Q3 = 2006, seuil = 0.5 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.72446 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.27554 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.2326555\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} +& \multicolumn{2}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 32779 & 11359 \\ + 1 & 2418 & 3444 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.931301\\ +& Taux de 'faux positifs' : 0.06869904\\ +& Taux de 'faux n\'egatifs' : 0.7673445\\ +\hline +\end{tabu}} + + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_test\_logic : Q1 = 1994, seuil = 0.5 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.80248 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.19752 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.935429\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} +& \multicolumn{2}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 3820 & 2506 \\ + 1 & 7370 & 36304 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.3413762\\ +& Taux de 'faux positifs' : 0.6586238\\ +& Taux de 'faux n\'egatifs' : 0.06457099\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_test\_logic : Me = 2002, seuil = 0.5 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.71346 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.28654 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7675815\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} +& \multicolumn{2}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 15732 & 6038 \\ + 1 & 8289 & 19941 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.6549269\\ +& Taux de 'faux positifs' : 0.3450731\\ +& Taux de 'faux n\'egatifs' : 0.2324185\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_test\_logic : Q3 = 2006, seuil = 0.5 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.71868 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.28132 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.2301975\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} +& \multicolumn{2}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 32496 & 11497 \\ + 1 & 2569 & 3438 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.9267361\\ +& Taux de 'faux positifs' : 0.07326394\\ +& Taux de 'faux n\'egatifs' : 0.7698025\\ +\hline +\end{tabu}} + +## Annexe 8 - Matrice de confusion pour la logistique multimodale avec seuils arbitraires {#a_8} +```{r, message=FALSE, warning=FALSE, include=FALSE} +analyze.perf.multi <- function(sample_categ, sample_categ_predict, print=1) { + + conf <- table(sample_categ_predict, sample_categ) + n <- length(conf[1,]) + + eff_lig <- array(dim = n) + eff_col <- array(dim = n) + eff_good_lig <- array(dim = n) + eff_good_col <- array(dim = n) + + for(i in 1:n) { + eff_lig[i] <- sum(conf[i,]) + eff_col[i] <- sum(conf[,i]) + eff_good_lig[i] <- conf[i,i]/eff_lig[i] + eff_good_col[i] <- conf[i,i]/eff_col[i] + } + acc <- sum(diag(conf))/sum(conf) ; err <- 1 - acc + + if(print==1) { + cat("\n") + cat("------------------------------------------------\n") + cat("Analyse des performances du modèles\n") + cat("Comparaison entre réponses prédites et attendues\n") + cat("------------------------------------------------\n") + cat("\n Matrice de confusion :\n") ; print(conf) + cat("\n Performance : ") ; cat(acc) + cat("\n Taux d'erreur : ") ; cat(err) + for(i in 1:n) { + cat("\n Taux BP a priori (mod ") ; cat(i) ; cat(") : ") ; cat(eff_good_col[i]) + cat(" | Taux BP a posteriori (mod ") ; cat(i) ; cat(") : ") ; cat(eff_good_lig[i]) + } + } + else{ + res <- list("conf" = conf, "acc" = acc, "err" = err, + "eff_lig" = eff_lig, "eff_col" = eff_col, + "eff_good_lig" = eff_good_lig, "eff_good_col" = eff_good_col) + return(res) + } +} + +sample_train_categ <- ifelse(sample_train$Year < q_fig[2], 1, + ifelse(sample_train$Year < q_fig[3], 2, + ifelse(sample_train$Year < q_fig[4], 3, 4))) + +sample_train_categ_predict <- ifelse(res_logic_Q2$fitted.values <= seuilQ2, + ifelse(res_logic_Q1$fitted.values <= seuilQ1, 1, 2), + ifelse(res_logic_Q3$fitted.values <= seuilQ3, 3, 4)) + +cat("sample_train_categ | seuils arbitraires \n") ; analyze.perf.multi(sample_train_categ,sample_train_categ_predict) + +sample_test_categ <- ifelse(sample_test$Year < q_fig[2], 1, + ifelse(sample_test$Year < q_fig[3], 2, + ifelse(sample_test$Year < q_fig[4], 3, 4))) + +predict_test_Q1 <- predict.glm(res_logic_Q1, sample_test_logic_Q1, type="response") +predict_test_Q2 <- predict.glm(res_logic_Q2, sample_test_logic_Q2, type="response") +predict_test_Q3 <- predict.glm(res_logic_Q3, sample_test_logic_Q3, type="response") + +sample_test_categ_predict <- ifelse(predict_test_Q2 <= seuilQ2, + ifelse(predict_test_Q1 <= seuilQ1, 1, 2), + ifelse(predict_test_Q3 <= seuilQ3, 3, 4)) + +cat("sample_test_categ | seuils arbitraires \n") ; analyze.perf.multi(sample_test_categ,sample_test_categ_predict) + +``` + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline +\textit{Analyse des performances du mod\`eles} & Taux BP a priori (mod 1) : 0.3602398 \tabularnewline +\textit{Comparaison entre r\'eponses pr\'edites et attendues} & Taux BP a priori (mod 2) : 0.3974267 \tabularnewline +\tabucline[on 2pt]{-1} +Matrice de confusion : & Taux BP a priori (mod 3) : 0.5805173 \\ +\tabulinesep=0.25mm{ +\multirow{4}{*}{ +\begin{tabu}{r|l|l|l|l} + & \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 1 & 2 & 3 & 4 \\ + 1 & 4146 & 1428 & 522 & 402 \\ + 2 & 5369 & 5004 & 2510 & 2380 \\ + 3 & 1861 & 5497 & 6442 & 8577 \\ + 4 & 133 & 662 & 1623 & 3444 +\end{tabu} +}} & Taux BP a priori (mod 4) : 0.2326555\\ +& Taux BP a posteriori (mod 1) : 0.6380425\\ +& Taux BP a posteriori (mod 2) : 0.3278517 \\ +& Taux BP a posteriori (mod 3) : 0.2878849\\ +\tabucline[on 2pt]{-1} +\textbf{Performance : 0.38072} & Taux BP a posteriori (mod 4) : 0.5875128\\ +\textbf{Taux d'erreur : 0.61928} & \\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline +\textit{Analyse des performances du mod\`eles} & Taux BP a priori (mod 1) : 0.3602398 \tabularnewline +\textit{Comparaison entre r\'eponses pr\'edites et attendues} & Taux BP a priori (mod 2) : 0.3974267 \tabularnewline +\tabucline[on 2pt]{-1} +Matrice de confusion : & Taux BP a priori (mod 3) : 0.5805173 \\ +\tabulinesep=0.25mm{ +\multirow{4}{*}{ +\begin{tabu}{r|l|l|l|l} + & \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 1 & 2 & 3 & 4 \\ + 1 & 4146 & 1428 & 522 & 402 \\ + 2 & 5369 & 5004 & 2510 & 2380 \\ + 3 & 1861 & 5497 & 6442 & 8577 \\ + 4 & 133 & 662 & 1623 & 3444 +\end{tabu} +}} & Taux BP a priori (mod 4) : 0.2326555\\ +& Taux BP a posteriori (mod 1) : 0.6380425\\ +& Taux BP a posteriori (mod 2) : 0.3278517 \\ +& Taux BP a posteriori (mod 3) : 0.2878849\\ +\tabucline[on 2pt]{-1} +\textbf{Performance : 0.38072} & Taux BP a posteriori (mod 4) : 0.5875128\\ +\textbf{Taux d'erreur : 0.61928} & \\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline +\textit{Analyse des performances du mod\`eles} &\large \bfseries sample\_test\_categ, seuils arbitraires \tabularnewline +\tabucline[on 2pt]{2} +\textit{Comparaison entre r\'eponses pr\'edites et attendues} & Taux BP a priori (mod 1) : 0.3412869 \tabularnewline +\tabucline[on 2pt]{-1} +Matrice de confusion : & Taux BP a priori (mod 2) : 0.3989557 \\ +\tabulinesep=0.25mm{ +\multirow{4}{*}{ +\begin{tabu}{r|l|l|l|l} + & \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 1 & 2 & 3 & 4 \\ + 1 & 3819 & 1411 & 590 & 504 \\ + 2 & 5383 & 5119 & 2421 & 2523 \\ + 3 & 1855 & 5616 & 6282 & 8470 \\ + 4 & 133 & 685 & 1751 & 3438 +\end{tabu} +}} & Taux BP a priori (mod 3) : 0.5688156\\ +& Taux BP a priori (mod 4) : 0.2301975\\ +& Taux BP a posteriori (mod 1) : 0.6038899 \\ +& Taux BP a posteriori (mod 2) : 0.3314127\\ +\tabucline[on 2pt]{-1} +\textbf{Performance : 0.37316} & Taux BP a posteriori (mod 3) : 0.2826801\\ +\textbf{Taux d'erreur : 0.62684} & Taux BP a posteriori (mod 4) : 0.5723323\\ +\hline +\end{tabu}} + +## Annexe 9 - Courbes ROC d'après la régression logistique {#a_9} +```{r Régression Logistique, echo=FALSE, fig.height=2, message=FALSE, warning=FALSE} +perf_roc <- function(res, tab, var_exp) { + pred <- prediction(predict.glm(res, tab, type="response"), as.numeric(t(tab[var_exp]))) + perf <- performance(pred,"tpr","fpr") + plot(perf,colorize=TRUE, asp = 1) + lines(c(0,1),c(0,1)) + return(perf) +} + +par(mfrow=c(1,3)) +perf_log_Q1 <- perf_roc(res_logic_Q1, sample_test_logic_Q1, "Year") +perf_log_Q2 <- perf_roc(res_logic_Q2, sample_test_logic_Q2, "Year") +perf_log_Q3 <- perf_roc(res_logic_Q3, sample_test_logic_Q3, "Year") +par(new=T) +``` + +## Annexe 10 - Matrice de confusion pour la logistique avec seuils pertinents {#a_10} +```{r, message=FALSE, warning=FALSE, include=FALSE} +search_cutoff_equal <- function(perf) { + n <- length(slot(perf,"x.values")[[1]]) + p <- floor(n/100) ; cat(p) + ind <- which((slot(perf,"y.values")[[1]][seq(1,n-p,p)] - slot(perf,"y.values")[[1]][seq(1+p,n,p)])/ + (slot(perf,"x.values")[[1]][seq(1,n-p,p)] - slot(perf,"x.values")[[1]][seq(1+p,n,p)])<1)[1] + cutoff <- slot(perf,"alpha.values")[[1]][(ind+1)*p] + return(cutoff) +} + +seuil_logic_Q1 <- search_cutoff_equal(perf_log_Q1) ; +seuil_logic_Q2 <- search_cutoff_equal(perf_log_Q2) ; +seuil_logic_Q3 <- search_cutoff_equal(perf_log_Q3) ; + +cat ("sample_train_logic : Q1 = ", q_fig[2], " | seuil = ", seuil_logic_Q1, "\n") ; analyze.perf(res_logic_Q1, "Year", seuil_logic_Q1) +cat ("sample_train_logic : Me = ", q_fig[3], " | seuil = ", seuil_logic_Q2, "\n") ; analyze.perf(res_logic_Q2, "Year", seuil_logic_Q2) +cat ("sample_train_logic : Q3 = ", q_fig[4], " | seuil = ", seuil_logic_Q3, "\n") ; analyze.perf(res_logic_Q3, "Year", seuil_logic_Q3) + +cat ("sample_test_logic : Q1 = ", q_fig[2], " | seuil = ", seuil_logic_Q1, "\n") ; analyze.perf(res_logic_Q1, "Year", seuil_logic_Q1, sample_test_logic_Q1) +cat ("sample_test_logic : Me = ", q_fig[3], " | seuil = ", seuil_logic_Q2, "\n") ; analyze.perf(res_logic_Q2, "Year", seuil_logic_Q2, sample_test_logic_Q2) +cat ("sample_test_logic : Q3 = ", q_fig[4], " | seuil = ", seuil_logic_Q3, "\n") ; analyze.perf(res_logic_Q3, "Year", seuil_logic_Q3, sample_test_logic_Q3) + +sample_train_categ <- ifelse(sample_train$Year < q_fig[2], 1, + ifelse(sample_train$Year < q_fig[3], 2, + ifelse(sample_train$Year < q_fig[4], 3, 4))) + +sample_train_categ_predict <- ifelse(res_logic_Q2$fitted.values <= seuil_logic_Q2, + ifelse(res_logic_Q1$fitted.values <= seuil_logic_Q1, 1, 2), + ifelse(res_logic_Q3$fitted.values <= seuil_logic_Q3, 3, 4)) + +cat ("sample_train_categ | seuils pertinents\n") ; analyze.perf.multi(sample_train_categ,sample_train_categ_predict) + +sample_test_categ <- ifelse(sample_test$Year < q_fig[2], 1, + ifelse(sample_test$Year < q_fig[3], 2, + ifelse(sample_test$Year < q_fig[4], 3, 4))) + +predict_test_Q1 <- predict.glm(res_logic_Q1, sample_test_logic_Q1, type="response") +predict_test_Q2 <- predict.glm(res_logic_Q2, sample_test_logic_Q2, type="response") +predict_test_Q3 <- predict.glm(res_logic_Q3, sample_test_logic_Q3, type="response") + +sample_test_categ_predict <- ifelse(predict_test_Q2 <= seuil_logic_Q2, + ifelse(predict_test_Q1 <= seuil_logic_Q1, 1, 2), + ifelse(predict_test_Q3 <= seuil_logic_Q3, 3, 4)) + +cat ("sample_test_categ | seuils pertinents\n") ; analyze.perf.multi(sample_test_categ,sample_test_categ_predict) + +``` + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_train\_logic : Q1 = 1994, seuil = 0.768732 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.76676 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.23324 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7663869\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} + \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 8839 & 8992 \\ + 1 & 2670 & 29499 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.7680076\\ +& Taux de 'faux positifs' : 0.2319924\\ +& Taux de 'faux n\'egatifs' : 0.2336131\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_train\_logic : Me = 2000, seuil = 0.5325237 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.72212 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.27788 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7342857\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} + \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 17088 & 6882 \\ + 1 & 7012 & 19018 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.7090456\\ +& Taux de 'faux positifs' : 0.2909544\\ +& Taux de 'faux n\'egatifs' : 0.2657143\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_train\_logic : Q3 = 2006, seuil = 0.2729759 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.64248 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.35752 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7779504\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} + \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 20608 & 3287 \\ + 1 & 14589 & 11516 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.5855044\\ +& Taux de 'faux positifs' : 0.4144956\\ +& Taux de 'faux n\'egatifs' : 0.2220496\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_test\_logic : Q3 = 1994, seuil = 0.768732 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.76478 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.23522 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7665035\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} + \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 8491 & 9062 \\ + 1 & 2699 & 29748 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.7588025\\ +& Taux de 'faux positifs' : 0.2411975\\ +& Taux de 'faux n\'egatifs' : 0.2334965\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_test\_logic : Me = 2002, seuil = 0.5325237 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.71368 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.28632 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7240848\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} + \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 16873 & 7168 \\ + 1 & 7148 & 18811 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.702427\\ +& Taux de 'faux positifs' : 0.297573\\ +& Taux de 'faux n\'egatifs' : 0.2759152\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline + \multicolumn{2}{|c|}{\large \bfseries sample\_test\_logic : Q3 = 2006, seuil = 0.2729759 } \tabularnewline +\tabucline[on 2pt]{1-2} +Analyse des performances du mod\`eles & Performance : 0.63522 \tabularnewline +Comparaison entre r\'eponses pr\'edites et attendues & Taux d'erreur : 0.36478 \tabularnewline +\tabucline[on 2pt]{1-2} +Matrice de confusion : & Sensibilit\'e (vrais positifs) : 0.7584868\tabularnewline +\multirow{2}{*}{ +\tabulinesep=0.5mm{ +\begin{tabu}{l|l|l} + \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 0 & 1 \\ + 0 & 20433 & 3607 \\ + 1 & 14632 & 11328 +\end{tabu} +}} & Sp\'ecificit\'e (vrais n\'egatifs) : 0.5827178\\ +& Taux de 'faux positifs' : 0.4172822\\ +& Taux de 'faux n\'egatifs' : 0.2415132\\ +\hline +\end{tabu}} + + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline +\textit{Analyse des performances du mod\`eles} &\large \bfseries sample\_train\_categ, seuils pertinents \tabularnewline +\tabucline[on 2pt]{2} +\textit{Comparaison entre r\'eponses pr\'edites et attendues} & Taux BP a priori (mod 1) : 0.7601877 \tabularnewline +\tabucline[on 2pt]{-1} +Matrice de confusion : & Taux BP a priori (mod 2) : 0.1983163 \\ +\tabulinesep=0.25mm{ +\multirow{4}{*}{ +\begin{tabu}{r|l|l|l|l} + & \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 1 & 2 & 3 & 4 \\ + 1 & 8749 & 4686 & 2088 & 1948 \\ + 2 & 1156 & 2497 & 1481 & 1365 \\ + 3 & 119 & 469 & 378 & 322 \\ + 4 & 1485 & 4939 & 7150 & 11168 +\end{tabu} +}} & Taux BP a priori (mod 3) : 0.03406326\\ +& Taux BP a priori (mod 4) : 0.7544417\\ +& Taux BP a posteriori (mod 1) : 0.5007727 \\ +& Taux BP a posteriori (mod 2) : 0.384213\\ +\tabucline[on 2pt]{-1} +\textbf{Performance : 0.45584} & Taux BP a posteriori (mod 3) : 0.2934783\\ +\textbf{Taux d'erreur : 0.54416} & Taux BP a posteriori (mod 4) : 0.4513782\\ +\hline +\end{tabu}} + +\tabulinesep=2mm{ +\begin{tabu}{|l|c|} +\hline +\textit{Analyse des performances du mod\`eles} &\large \bfseries sample\_test\_categ, seuils pertinents \tabularnewline +\tabucline[on 2pt]{2} +\textit{Comparaison entre r\'eponses pr\'edites et attendues} & Taux BP a priori (mod 1) : 0.7522788 \tabularnewline +\tabucline[on 2pt]{-1} +Matrice de confusion : & Taux BP a priori (mod 2) : 0.2058296 \\ +\tabulinesep=0.25mm{ +\multirow{4}{*}{ +\begin{tabu}{r|l|l|l|l} + & \multicolumn{3}{r}{Valeur de r\'ef\'erence} \\ + Valeur pr\'edite & 1 & 2 & 3 & 4 \\ + 1 & 8418 & 4637 & 2095 & 2097 \\ + 2 & 178 & 2641 & 1420 & 1556 \\ + 3 & 156 & 490 & 396 & 354 \\ + 4 & 1438 & 5063 & 7133 & 10928 +\end{tabu} +}} & Taux BP a priori (mod 3) : 0.03585657\\ +& Taux BP a priori (mod 4) : 0.7317041\\ +& Taux BP a posteriori (mod 1) : 0.4880849 \\ +& Taux BP a posteriori (mod 2) : 0.3886681\\ +\tabucline[on 2pt]{-1} +\textbf{Performance : 0.44766} & Taux BP a posteriori (mod 3) : 0.2836676\\ +\textbf{Taux d'erreur : 0.55234} & Taux BP a posteriori (mod 4) : 0.4449149\\ +\hline +\end{tabu}} + +## Annexe 11 - Représentation de l'ACP {#a_11} +```{r PCA, echo=FALSE, fig.height=2, message=FALSE, warning=FALSE} +library("FactoMineR") +library("factoextra") + +<<<<<<< HEAD +======= + +>>>>>>> MRR presque finie +res.pca <- PCA(sample_train[,-1], graph = FALSE) + +# get_eigenvalue(res.pca) +fviz_eig(res.pca, addlabels = TRUE) + +ind <- get_pca_ind(res.pca) ; #ind +var <- get_pca_var(res.pca) ; #var +fviz_pca_ind(res.pca, + geom.ind = c("point"), + col.ind = sample_train$Year, # colorer by groups + legend.title = "Groups") + +fviz_pca_var(res.pca, geom.var = c("arrow"), col.var = "contrib", + gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")) + +# 20 % expliqué dans le plan => l'impression que les variables sont fortement non corrélées +``` \ No newline at end of file diff --git a/Rapport Final/Rapport_Final.pdf b/Rapport Final/Rapport_Final.pdf new file mode 100755 index 0000000000000000000000000000000000000000..af6aa8149040a057bd6b0ec46ccbc3333ddee578 Binary files /dev/null and b/Rapport Final/Rapport_Final.pdf differ diff --git a/Soutenance/.DS_Store b/Soutenance/.DS_Store new file mode 100755 index 0000000000000000000000000000000000000000..56d40585678ca03cb1007f3473a722bbe390c792 Binary files /dev/null and b/Soutenance/.DS_Store differ diff --git a/Soutenance/Soutenance.pdf b/Soutenance/Soutenance.pdf new file mode 100755 index 0000000000000000000000000000000000000000..a524f3a7e33bb4ad3411c53e1ee3a8856d8a46a2 Binary files /dev/null and b/Soutenance/Soutenance.pdf differ diff --git "a/Soutenance/images/ConclusionMod\303\250lesLin\303\251aires.png" "b/Soutenance/images/ConclusionMod\303\250lesLin\303\251aires.png" new file mode 100755 index 0000000000000000000000000000000000000000..2969668166dbba7633b2a90bad9ecfdf47e0fe0c Binary files /dev/null and "b/Soutenance/images/ConclusionMod\303\250lesLin\303\251aires.png" differ diff --git "a/Soutenance/images/Distribution des chansons selon leur ann\303\251e de sortie.png" "b/Soutenance/images/Distribution des chansons selon leur ann\303\251e de sortie.png" new file mode 100755 index 0000000000000000000000000000000000000000..9ca6971a8ab76276243c63ee2bf1b48dfbe1b415 Binary files /dev/null and "b/Soutenance/images/Distribution des chansons selon leur ann\303\251e de sortie.png" differ diff --git a/Soutenance/images/PCA.png b/Soutenance/images/PCA.png new file mode 100755 index 0000000000000000000000000000000000000000..3a93718fe56e6dd9139b7d1311a0a207bef5f3d9 Binary files /dev/null and b/Soutenance/images/PCA.png differ diff --git a/Soutenance/images/PCA2.png b/Soutenance/images/PCA2.png new file mode 100755 index 0000000000000000000000000000000000000000..bc8c5f17125d4a801fbe77fb13f8e268b1275fdd Binary files /dev/null and b/Soutenance/images/PCA2.png differ diff --git a/Soutenance/images/ROC.png b/Soutenance/images/ROC.png new file mode 100755 index 0000000000000000000000000000000000000000..eec2abd25bd0f130814b9e752042f5a8c494c825 Binary files /dev/null and b/Soutenance/images/ROC.png differ diff --git a/Soutenance/images/Rplot.png b/Soutenance/images/Rplot.png new file mode 100755 index 0000000000000000000000000000000000000000..1cec036ca74187ad55d60a1cf48efeb044f31a55 Binary files /dev/null and b/Soutenance/images/Rplot.png differ diff --git a/YearPredictionMSD.txt b/YearPredictionMSD.txt new file mode 100755 index 0000000000000000000000000000000000000000..55c9aac015296d47fc92d01648b11e04282794b5 Binary files /dev/null and b/YearPredictionMSD.txt differ diff --git a/code.R b/code.R new file mode 100755 index 0000000000000000000000000000000000000000..6436ce015183ad6995c45789647b4084f93157cf --- /dev/null +++ b/code.R @@ -0,0 +1,737 @@ +####################################################################################### +####################################################################################### + +# Projet - YearPredictionMSD Data Set + +# Source : +# Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. +# Irvine, CA: University of California, School of Information and Computer Science. + +# Code : +# Thibault CORDIER & Julien KHAMPHOUSONE + +####################################################################################### +####################################################################################### + +setwd("~/Documents/Projet/Projets-ENSIIE-MA/MRR") +rm(list=ls()) +graphics.off() + +require(ggplot2) +require(GGally) + +#lm.ridge +require(MASS) + +#glmnet +require(glmnet) +require(corrplot) + +#roc +require(ROCR) + + + + + +####################################################################################### +# Initialisation de notre base de données +####################################################################################### + +####################################################################################### +# Renommages des variables de nos données +my_colnames <- rep(NA, 91) +my_colnames[1] <- "Year" +for(i in (1:12)) my_colnames[i+1] <- paste("T_Ave_",i,sep="") +for(i in (1:78)) my_colnames[i+13] <- paste("T_Cov_",i,sep="") + +# Importation de nos données d'entrainement et de test +t0 <- proc.time() +data_train <- read.table("YearPredictionMSD.txt", sep=",", nrows = 463715, col.names = my_colnames) +tf <- proc.time() ; tf-t0 ; t0 <- proc.time() +data_test <- read.table("YearPredictionMSD.txt", sep=",", col.names = my_colnames, skip = 463715) +tf <- proc.time() ; tf-t0 + +sample_train <- data_train[sample(n_train,50000),] +sample_test <- data_test[sample(n_test,50000),] + +####################################################################################### +# Indicateurs +sum(is.na(data_train[,1])) + + + + + +######################################################################################## +# Visualisation des données +######################################################################################## + +# ggpairs(sample_train[,c(1,sample(2:91,5))]) + +# ggplot(sample_train, aes(Year, T_Ave_1)) + geom_point() + ggtitle("Scatterplot : +# Year en fonction de T_Ave_1") + +# { i <- floor(runif(1, 2, 91)) ; +# ggplot(sample_train, aes_string("Year", my_colnames[i])) + geom_point() + +# ggtitle(paste("Scatterplot : Year = f(", my_colnames[i], ")", sep="")) +# ggally_box_no_facet(sample_train, aes_string("Year", my_colnames[i])) +# } + +summary(data_train["Year"]) + +ggplot(data_train, aes(Year)) + geom_bar() + + ggtitle("Barplot : Distribution des chansons selon leur année de sortie") + +ggplot(data_train, aes(1,Year)) + geom_boxplot() + + ggtitle("Boxplot : Distribution des chansons selon leur année de sortie") + + + + + +######################################################################################## +# Construction d’un premier modèle de régression linéaire +######################################################################################## + +reg <- lm(Year~., sample_train) +plot(reg$model$Year, reg$residuals, xlab = "Year", ylab = "Résidu") + +mse_reg <- sqrt(mean((reg$model$Year-reg$fitted.values)^2)) ; mse_reg + + + + + +######################################################################################## +# Discussion sur une filtration de nos données d’entrainement +######################################################################################## +# Sélection partielle de nos données d'entrainement et de test +n_train <- dim(data_train)[1] +n_test <- dim(data_test)[1] + +# Sélection au préalable des observations (outliers) +1994 - 1.5*(2006-1994) # Méthode des moustaches : 1976 +# 1970 -> 2010 + +selec_obs <- function(db,n) { + new_db <- NULL + for(y in 1970:2010) { + indk <- which(db[,"Year"]==y) + new_db <- rbind(new_db, db[sample(indk,min(n,length(indk))),]) + } + new_db +} + +sample_train_by_year <- selec_obs(data_train,1000) # sample_train contient 41000 observations +sample_test_by_year <- selec_obs(data_test,1000) + +ggplot(sample_train, aes(Year)) + geom_bar() + + ggtitle("Barplot : Distribution des chansons selon leur année de sortie après filtration") + + + + + +######################################################################################## +# Modèle de régression linéaire : Analyse des performances selon la filtration des observations +######################################################################################## +# Modèles +res_all <- glm(Year~., data=sample_train) +res_all_by_year <- glm(Year~., data=sample_train_by_year) + +# Prédictions +predict_all <- predict.glm(res_all, newdata = sample_test) +predict_all_by_year <- predict.glm(res_all_by_year,newdata = sample_test) +predict_all_2 <- predict.glm(res_all, newdata = sample_test_by_year) +predict_all_by_year_2 <- predict.glm(res_all_by_year,newdata = sample_test_by_year) + +# Résidus +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(res_all$model$Year, res_all$residuals, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Train", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(res_all_by_year$model$Year, res_all_by_year$residuals, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Train", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +plot(sample_test$Year, sample_test$Year - predict_all, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Test1", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(sample_test$Year, sample_test$Year - predict_all_by_year, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Test1", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +plot(sample_test_by_year$Year, sample_test_by_year$Year - predict_all_2, + xlab = "Y", ylab = "résidus", main="Résidus ; Sans sélection ; Test2", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(sample_test_by_year$Year, sample_test_by_year$Year - predict_all_by_year_2, + xlab = "Y", ylab = "résidus by year", main="Résidus ; Avec sélection ; Test2", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +# Y par rapport à Y^ +par(mfrow=c(2,2),mar = rep(2, 4)) +plot(sample_test$Year, predict_all, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Sans sélection ; Test1", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test$Year, predict_all_by_year, + xlab = "Y", ylab = "Y^ by year", main="Y/Y^ ; Avec sélection ; Test1", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test_by_year$Year, predict_all_2, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Sans sélection ; Test2", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +plot(sample_test_by_year$Year, predict_all_by_year_2, + xlab = "Y", ylab = "Y^ by year", main="Y/Y^ ; Avec sélection ; Test2", cex = 0.3, + xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) + +# MSE +mse_all <- sqrt(mean((res_all$model$Year - res_all$fitted.values)^2)) ; mse_all +mse_year <- sqrt(mean((res_all_by_year$model$Year - res_all_by_year$fitted.values)^2)) ; mse_year + +mse_all_p1 <- sqrt(mean((predict_all - sample_test$Year)^2)) ; mse_all_p1 +mse_year_p1 <- sqrt(mean((predict_all_by_year - sample_test$Year)^2)) ; mse_year_p1 + +mse_all_p2 <- sqrt(mean((predict_all_2 - sample_test_by_year$Year)^2)) ; mse_all_p2 +mse_year_p2 <- sqrt(mean((predict_all_by_year_2 - sample_test_by_year$Year)^2)) ; mse_year_p2 + +# KFold +kfold <- function(tab,var_exp,k) { + n <- nrow(tab) + m <- n %/% k + X <- tab[sample(nrow(tab)),] + error_test <- c() + error_train <- c() + E_train <- 0 + E_test <- 0 + + for(i in 1:k){ + ind_sample = ((i-1)*m+1):(i*m) + mod = glm(as.formula(paste(var_exp,"~.")), data=X[-ind_sample,]) + + pY = predict.glm(mod, X[-ind_sample,]) + E_train = E_train + sqrt(mean((X[-ind_sample,var_exp] - pY)^2)) + error_train <- c(error_train,sqrt(mean((X[-ind_sample,var_exp] - pY)^2))) + + pY = predict.glm(mod, X[ind_sample,]) + E_test = E_test + sqrt(mean((X[ind_sample,var_exp] - pY)^2)) + error_test <- c(error_test,sqrt(mean((X[ind_sample,var_exp] - pY)^2))) + } + E_train <- E_train/k + E_test <- E_test/k + + res <- list("Etrain" = E_train, "Etest" = E_test, + "ErrorTrain" = error_train, "ErrorTest" = error_test, + "n" = n, "m" = m, "k" = k) + return(res) +} + +mse_kf <- kfold(sample_train, "Year", 10) ; mse_kf +mse_kf_by_year <- kfold(sample_train_by_year, "Year", 10) ; mse_kf_by_year + + + + + +######################################################################################## +# Modèle de régression linéaire avec sélection de variables “forward” +######################################################################################## + +# res_all <- glm(Year~., data=sample_train) +# res_nul <- glm(Year~1, data=sample_train) +# fw_all <- step(res_nul, list(upper=res_all), direction="forward") +# +# predict_all <- predict.glm(res_all, newdata = sample_test) +# predict_all_fw <- predict.glm(fw_all, newdata = sample_test) +# +# # Résidus +# par(mfrow=c(2,2),mar = rep(2, 4)) +# plot(res_all$model$Year, res_all$residuals, +# xlab = "Y", ylab = "résidus", main="Résidus ; Standard ; Train", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# plot(fw_all$model$Year, fw_all$residuals, +# xlab = "Y", ylab = "résidus by year", main="Résidus ; Forward ; Train", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# +# plot(sample_test$Year, sample_test$Year - predict_all, +# xlab = "Y", ylab = "résidus", main="Résidus ; Standard ; Test", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# plot(sample_test$Year, sample_test$Year - predict_all_fw, +# xlab = "Y", ylab = "résidus", main="Résidus ; Forward ; Test", cex = 0.3, +# xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +# +# # Y par rapport à Y^ +# par(mfrow=c(1,2),mar = rep(2, 4)) +# plot(sample_test$Year, predict_all, +# xlab = "Y", ylab = "Y^", main="Y/Y^ ; Standard ; Test", cex = 0.3, +# xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +# plot(sample_test$Year, predict_all_fw, +# xlab = "Y", ylab = "Y^", main="Y/Y^ ; Forward ; Test", cex = 0.3, +# xlim=c(1920, 2030), ylim=c(1960, 2030)) ; abline(a=0,b=1) +# +# # MSE +# sqrt(mean((res_all$model$Year - res_all$fitted.values)^2)) # 9.624977 +# sqrt(mean((fw_all$model$Year - fw_all$fitted.values)^2)) # 9.626477 +# +# sqrt(mean((predict_all - sample_test$Year)^2)) # 9.45553 +# sqrt(mean((predict_all_fw - sample_test$Year)^2)) # 9.457589 + + + + + +####################################################################################### +# Recherche de variables impactants les résidus +####################################################################################### + +####################################################################################### +# Corrélation +require(corrplot) +graphics.off() +par(mfrow = c(1,1)) +corrplot(cor(sample_train)[1:30,1:30]) + +plot(res_all$model$T_Ave_6, res_all$residuals, xlab = "Year", ylab = "Résidu") + + + + + +######################################################################################## +# Modèle de régression régularisée Lasso-Ridge-Elastic +######################################################################################## + +x.train<- as.matrix(scale(sample_train[2:length(sample_train)])) +x.test <- as.matrix(scale(sample_test[2:length(sample_test)])) + +y.train <- sample_train[,c("Year")] +y.test <- sample_test[,c("Year")] + +for (i in 0:4) { + assign(paste("fit", i, sep=""), cv.glmnet(x.train, y.train, type.measure="mse", + alpha=i/4, family="gaussian")) + print(i); +} + +#fit0$lambda.min ; fit1$lambda.min ; fit2$lambda.min ; fit3$lambda.min ; fit4$lambda.min + +#plot(fit4, main="LASSO") +#plot(fit3, main="Elastic Net 0.75") +#plot(fit2, main="Elastic Net 0.5") +#plot(fit1, main="Elastic Net 0.25") +#plot(fit0, main="Ridge") + +# Modèles +res_all <- glm(Year~., data=sample_train) +fit.ridge <- glmnet(x.train, y.train, family = "gaussian", alpha=0, lambda = fit0$lambda.min) +fit.elnet0.25 <- glmnet(x.train, y.train, family = "gaussian", alpha=0.25, lambda = fit1$lambda.min) +fit.elnet0.5 <- glmnet(x.train, y.train, family = "gaussian", alpha=0.5, lambda = fit2$lambda.min) +fit.elnet0.75 <- glmnet(x.train, y.train, family = "gaussian", alpha=0.75, lambda = fit3$lambda.min) +fit.lasso <- glmnet(x.train, y.train, family = "gaussian", alpha=1, lambda = fit4$lambda.min) + +# Prédictions +yhatall <- predict.glm(res_all, newdata = sample_test) +yhat0 <- predict(fit.ridge, s=fit0$lambda.min, newx=x.test) +yhat1 <- predict(fit.elnet0.25, s=fit1$lambda.min, newx=x.test) +yhat2 <- predict(fit.elnet0.5, s=fit2$lambda.min, newx=x.test) +yhat3 <- predict(fit.elnet0.75, s=fit3$lambda.min, newx=x.test) +yhat4 <- predict(fit.lasso, s=fit4$lambda.min, newx=x.test) + +# Résidus +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(y.test, y.test - yhatall, + xlab = "Y", ylab = "résidus", main="Résidus ; Standard ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat0, + xlab = "Y", ylab = "résidus", main="Résidus ; Ridge ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat1, + xlab = "Y", ylab = "résidus", main="Résidus ; ElNet.0,25 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat2, + xlab = "Y", ylab = "résidus", main="Résidus ; ElNet.0,50 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat3, + xlab = "Y", ylab = "résidus", main="Résidus ; ElNet.0,75 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) +plot(y.test, y.test - yhat4, + xlab = "Y", ylab = "résidus", main="Résidus ; Lasso ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(-60, 40)) ; abline(a=0,b=0) + +# Y par rapport à Y^ +par(mfrow=c(3,2),mar = rep(2, 4)) +plot(y.test, yhatall, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Standard ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat0, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Ridge ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat1, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; ElNet.0,25 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat2, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; ElNet.0,50 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat3, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; ElNet.0,75 ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) +plot(y.test, yhat4, + xlab = "Y", ylab = "Y^", main="Y/Y^ ; Lasso ; Test", cex = 0.3, + xlim=c(1920, 2011), ylim=c(1960, 2031)) ; abline(a=0,b=1) + +# MSE +mse_all <- sqrt(mean((y.test-yhatall)^2)) +mse_ridge <- sqrt(mean((y.test-yhat0)^2)) +mse_elastic0.25 <- sqrt(mean((y.test-yhat1)^2)) +mse_elastic0.5 <- sqrt(mean((y.test-yhat2)^2)) +mse_elastic0.75 <- sqrt(mean((y.test-yhat3)^2)) +mse_lasso <- sqrt(mean((y.test-yhat4)^2)) + +mse_all ; mse_ridge ; mse_elastic0.25 ; mse_elastic0.5 ; mse_elastic0.75 ; mse_lasso + + + + + +######################################################################################## +# Prise en compte de la marge d’erreur sur la performance des modèles +######################################################################################## + +# # install.packages("gridExtra") +# require(gridExtra) +# grid.arrange(q1,q2,q1,q2) + +Yhat <- res_all$fitted.values +Yhat_test <- predict.glm(res_all, newdata = sample_test) +Y <- res_all$model$Year + +calcul_delta <- function(Y, Yhat, delta) { + n <- length(Y) + res <- 0 + for(i in 1:n) if(abs(Yhat[i] - Y[i]) < delta) res = res + 1 + res/n +} + +calcul_delta(Y,Yhat,20) +calcul_delta(Y,Yhat,15) +calcul_delta(Y,Yhat,10) +calcul_delta(Y,Yhat,5) +calcul_delta(Y,Yhat,2) + +V_delta <- seq(0,20,0.1) +V_Suc <- sapply(V_delta, calcul_delta, Y=Y, Yhat=Yhat) +q1 <- qplot(V_delta, V_Suc) + +calcul_delta(Y,Yhat_test,20) +calcul_delta(Y,Yhat_test,15) +calcul_delta(Y,Yhat_test,10) +calcul_delta(Y,Yhat_test,5) +calcul_delta(Y,Yhat_test,2) + +V_Suc_test <- sapply(V_delta, calcul_delta, Y=Y, Yhat=Yhat_test) +q2 <- qplot(V_delta, V_Suc_test) + +data_marge <- data.frame(cbind(c(V_delta,V_delta), c(V_Suc,V_Suc_test), + c(rep(1,length(V_delta)),rep(2,length(V_delta))))) +names(data_marge) <- c("Marge_Erreur","Bonnes_Predictions","Categorie") +qplot(data_marge$Marge_Erreur,data_marge$Bonnes_Predictions, color = data_marge$Categorie, + xlab = "Marge d'erreur en années", ylab = "Pourcentage de bonnes prédictions") + + + + + +######################################################################################## +# PCA +######################################################################################## +# install.packages(c("FactoMineR", "factoextra")) +library("FactoMineR") +library("factoextra") + +graphics.off() +par(mfrow=c(1,2)) +res.pca <- PCA(sample_train[,-1], graph = FALSE) + +# get_eigenvalue(res.pca) +fviz_eig(res.pca, addlabels = TRUE) + +ind <- get_pca_ind(res.pca) ; #ind +var <- get_pca_var(res.pca) ; #var + +fviz_pca_ind(res.pca, + geom.ind = c("point"), + col.ind = sample_train$Year, # colorer by groups + legend.title = "Groups") +fviz_pca_var(res.pca, geom.var = c("arrow"), col.var = "contrib", + gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")) + +# 20 % expliqué dans le plan => l'impression que les variables sont fortement non corrélées + + + + + +######################################################################################## +# Régression Logistique +######################################################################################## + +######################################################################################## +# Initialisation des données pour la logistique + +sample_train <- data_train[sample(n_train,50000),] +sample_test <- data_test[sample(n_test,50000),] + +q_fig <- quantile(data_train$Year) + +sample_train_logic_Q1 <- sample_train +sample_train_logic_Q2 <- sample_train +sample_train_logic_Q3 <- sample_train + +sample_train_logic_Q1$Year <- ifelse(sample_train_logic_Q1$Year >= q_fig[2], 1, 0) +sample_train_logic_Q2$Year <- ifelse(sample_train_logic_Q2$Year >= q_fig[3], 1, 0) +sample_train_logic_Q3$Year <- ifelse(sample_train_logic_Q3$Year >= q_fig[4], 1, 0) + +sample_test_logic_Q1 <- sample_test +sample_test_logic_Q2 <- sample_test +sample_test_logic_Q3 <- sample_test + +sample_test_logic_Q1$Year <- ifelse(sample_test_logic_Q1$Year >= q_fig[2], 1, 0) +sample_test_logic_Q2$Year <- ifelse(sample_test_logic_Q2$Year >= q_fig[3], 1, 0) +sample_test_logic_Q3$Year <- ifelse(sample_test_logic_Q3$Year >= q_fig[4], 1, 0) + +######################################################################################## +# Construction des modèles logistiques + +res_logic_Q1 <- glm(Year~., data=sample_train_logic_Q1, family=binomial) +res_logic_Q2 <- glm(Year~., data=sample_train_logic_Q2, family=binomial) +res_logic_Q3 <- glm(Year~., data=sample_train_logic_Q3, family=binomial) + +######################################################################################## +# La fonction analysze.perf affiche les éléments liés à la matrice de confusion +analyze.perf <- function(res, var, seuil, data = res$model, print = 1) { + + val_ref <- t(data[var]) + val_pred <- ifelse(predict.glm(res, newdata = data, type="response")>= seuil, 1, 0) + + conf <- table(val_pred, val_ref) + + TN <- conf[1,1] ; FP <- conf[2,1] ; FN <- conf[1,2] ; TP <- conf[2,2] + av_N <- FP + TN ; av_P <- FN + TP ; ap_N <- TN + FN ; ap_P <- FP + TP + n <- TN + FP + FN + TP + + vp <- TP / av_P ; fp <- FP / av_N + vn <- TN / av_N ; fn <- FN / av_P + err <- (FN + FP) / n ; acc <- (TP + TN) / n + + if(print==1) { + cat("\n") + cat("------------------------------------------------\n") + cat("Analyse des performances du modèles\n") + cat("Comparaison entre réponses prédites et attendues\n") + cat("------------------------------------------------\n") + cat("\n Matrice de confusion :\n") ; print(conf) + cat("\n Performance : ") ; cat(acc) + cat("\n Taux d'erreur : ") ; cat(err) + cat("\n Sensibilité (vrais positifs) : ") ; cat(vp) + cat("\n Spécificité (vrais négatifs) : ") ; cat(vn) + cat("\n Taux de 'faux positifs' : ") ; cat(fp) + cat("\n Taux de 'faux négatifs' : ") ; cat(fn) + } + else{ + res <- list("conf" = conf, "acc" = acc, "err" = err, + "vp" = vp, "vn" = vn, "fp" = fp, "fn" = fn) + return(res) + } +} + +######################################################################################## +# La fonction analysze.perf.mult affiche les éléments liés à la matrice de confusion +# pour plusieurs modalités +analyze.perf.multi <- function(sample_categ, sample_categ_predict, print=1) { + + conf <- table(sample_categ_predict, sample_categ) + n <- length(conf[1,]) + + eff_lig <- array(dim = n) + eff_col <- array(dim = n) + eff_good_lig <- array(dim = n) + eff_good_col <- array(dim = n) + + for(i in 1:n) { + eff_lig[i] <- sum(conf[i,]) + eff_col[i] <- sum(conf[,i]) + eff_good_lig[i] <- conf[i,i]/eff_lig[i] + eff_good_col[i] <- conf[i,i]/eff_col[i] + } + acc <- sum(diag(conf))/sum(conf) ; err <- 1 - acc + + if(print==1) { + cat("\n") + cat("------------------------------------------------\n") + cat("Analyse des performances du modèles\n") + cat("Comparaison entre réponses prédites et attendues\n") + cat("------------------------------------------------\n") + cat("\n Matrice de confusion :\n") ; print(conf) + cat("\n Performance : ") ; cat(acc) + cat("\n Taux d'erreur : ") ; cat(err) + for(i in 1:n) { + cat("\n Taux BP a priori (mod ") ; cat(i) ; cat(") : ") ; cat(eff_good_col[i]) + cat(" | Taux BP a posteriori (mod ") ; cat(i) ; cat(") : ") ; cat(eff_good_lig[i]) + } + } + else{ + res <- list("conf" = conf, "acc" = acc, "err" = err, + "eff_lig" = eff_lig, "eff_col" = eff_col, + "eff_good_lig" = eff_good_lig, "eff_good_col" = eff_good_col) + return(res) + } +} + +######################################################################################## +# Définition des seuils pour la prédiction +seuil <- 0.5 +seuilQ1 <- 0.5 +seuilQ2 <- 0.5 +seuilQ3 <- 0.5 + +######################################################################################## +# Réalisation des prédictions selon un seuil arbitraire : +# Modèle bimodal sur données d'entrainement +cat ("sample_train_logic : Q1 = ") ; cat(q_fig[2]) ; cat (" | seuil = ") ; cat(seuilQ1) ; cat("\n") ; analyze.perf(res_logic_Q1, "Year", seuilQ1) +cat ("sample_train_logic : Me = ") ; cat(q_fig[3]) ; cat (" | seuil = ") ; cat(seuilQ2) ; cat("\n") ; analyze.perf(res_logic_Q2, "Year", seuilQ2) +cat ("sample_train_logic : Q3 = ") ; cat(q_fig[4]) ; cat (" | seuil = ") ; cat(seuilQ3) ; cat("\n") ;analyze.perf(res_logic_Q3, "Year", seuilQ3) + +# Réalisation des prédictions selon un seuil arbitraire +# Modèle bimodal sur données de test +cat ("sample_test_logic : Q1 = ") ; cat(q_fig[2]) ; cat (" | seuil = ") ; cat(seuilQ1) ; cat("\n") ; analyze.perf(res_logic_Q1, "Year", seuilQ1, sample_test_logic_Q1) +cat ("sample_test_logic : Me = ") ; cat(q_fig[3]) ; cat (" | seuil = ") ; cat(seuilQ2) ; cat("\n") ; analyze.perf(res_logic_Q2, "Year", seuilQ2, sample_test_logic_Q2) +cat ("sample_test_logic : Q3 = ") ; cat(q_fig[4]) ; cat (" | seuil = ") ; cat(seuilQ3) ; cat("\n") ; analyze.perf(res_logic_Q3, "Year", seuilQ3, sample_test_logic_Q3) + +# APQ1 <- analyze.perf(res_logic_Q1, "Year", seuilQ1, print=FALSE) +# APQ2 <- analyze.perf(res_logic_Q2, "Year", seuilQ2, print=FALSE) +# APQ3 <- analyze.perf(res_logic_Q3, "Year", seuilQ3, print=FALSE) + +######################################################################################## +# Réalisation des prédictions selon un seuil arbitraire +# Modèle multimodal sur données d'entrainement + +sample_train_categ <- ifelse(sample_train$Year < q_fig[2], 1, + ifelse(sample_train$Year < q_fig[3], 2, + ifelse(sample_train$Year < q_fig[4], 3, 4))) + +sample_train_categ_predict <- ifelse(res_logic_Q2$fitted.values <= seuilQ2, + ifelse(res_logic_Q1$fitted.values <= seuilQ1, 1, 2), + ifelse(res_logic_Q3$fitted.values <= seuilQ3, 3, 4)) + +cat ("sample_train_categ | seuils arbitraires") ; cat("\n") ; analyze.perf.multi(sample_train_categ,sample_train_categ_predict) + +# Réalisation des prédictions selon un seuil arbitraire +# Modèle multimodal sur données de test + +sample_test_categ <- ifelse(sample_test$Year < q_fig[2], 1, + ifelse(sample_test$Year < q_fig[3], 2, + ifelse(sample_test$Year < q_fig[4], 3, 4))) + +predict_test_Q1 <- predict.glm(res_logic_Q1, sample_test_logic_Q1, type="response") +predict_test_Q2 <- predict.glm(res_logic_Q2, sample_test_logic_Q2, type="response") +predict_test_Q3 <- predict.glm(res_logic_Q3, sample_test_logic_Q3, type="response") + +sample_test_categ_predict <- ifelse(predict_test_Q2 <= seuilQ2, + ifelse(predict_test_Q1 <= seuilQ1, 1, 2), + ifelse(predict_test_Q3 <= seuilQ3, 3, 4)) + +cat ("sample_test_categ | seuils arbitraires") ; cat("\n") ; analyze.perf.multi(sample_test_categ,sample_test_categ_predict) + +######################################################################################## +# Visulaisation de la courbe ROC + +perf_roc <- function(res, tab, var_exp) { + pred <- prediction(predict.glm(res, tab, type="response"), as.numeric(t(tab[var_exp]))) + perf <- performance(pred,"tpr","fpr") + plot(perf,colorize=TRUE, asp = 1) + lines(c(0,1),c(0,1)) + return(perf) +} + +graphics.off() +par(mfrow=c(1,3)) +perf_log_Q1 <- perf_roc(res_logic_Q1, sample_test_logic_Q1, "Year") +perf_log_Q2 <- perf_roc(res_logic_Q2, sample_test_logic_Q2, "Year") +perf_log_Q3 <- perf_roc(res_logic_Q3, sample_test_logic_Q3, "Year") +par(new=T) + +######################################################################################## +# Calcul de l'AUC + +auc_roc <- function(res, tab, var_exp) { + pred <- prediction(predict.glm(res, tab, type="response"), as.numeric(t(tab[var_exp]))) + perf <- performance(pred,"tpr","fpr", measure = "auc") + return(slot(perf, "y.values")[[1]]) +} + +perf_auc_Q1 <- auc_roc(res_logic_Q1, sample_test_logic_Q1 , "Year") +perf_auc_Q1 +perf_auc_Q2 <- auc_roc(res_logic_Q2, sample_test_logic_Q2, "Year") +perf_auc_Q2 +perf_auc_Q3 <- auc_roc(res_logic_Q3, sample_test_logic_Q3, "Year") +perf_auc_Q3 + +######################################################################################## +# Recherche de seuils pertinents + +search_cutoff_equal <- function(perf) { + n <- length(slot(perf,"x.values")[[1]]) + p <- floor(n/100) ; cat(p) + ind <- which((slot(perf,"y.values")[[1]][seq(1,n-p,p)] - slot(perf,"y.values")[[1]][seq(1+p,n,p)])/ + (slot(perf,"x.values")[[1]][seq(1,n-p,p)] - slot(perf,"x.values")[[1]][seq(1+p,n,p)])<1)[1] + cutoff <- slot(perf,"alpha.values")[[1]][(ind+1)*p] + return(cutoff) +} + +seuil_logic_Q1 <- search_cutoff_equal(perf_log_Q1) ; seuil_logic_Q1 +seuil_logic_Q2 <- search_cutoff_equal(perf_log_Q2) ; seuil_logic_Q2 +seuil_logic_Q3 <- search_cutoff_equal(perf_log_Q3) ; seuil_logic_Q3 + +######################################################################################## +# Réalisation des prédictions selon un seuil pertinent : +# Modèle bimodal sur données d'entrainement +cat ("sample_train_logic : Q1 = ") ; cat(q_fig[2]) ; cat (" | seuil = ") ; cat(seuil_logic_Q1) ; cat("\n") ; analyze.perf(res_logic_Q1, "Year", seuil_logic_Q1) +cat ("sample_train_logic : Me = ") ; cat(q_fig[3]) ; cat (" | seuil = ") ; cat(seuil_logic_Q2) ; cat("\n") ; analyze.perf(res_logic_Q2, "Year", seuil_logic_Q2) +cat ("sample_train_logic : Q3 = ") ; cat(q_fig[4]) ; cat (" | seuil = ") ; cat(seuil_logic_Q3) ; cat("\n") ; analyze.perf(res_logic_Q3, "Year", seuil_logic_Q3) + +# Réalisation des prédictions selon un seuil pertinent : +# Modèle bimodal sur données de test +cat ("sample_test_logic : Q1 = ") ; cat(q_fig[2]) ; cat (" | seuil = ") ; cat(seuil_logic_Q1) ; cat("\n") ; analyze.perf(res_logic_Q1, "Year", seuil_logic_Q1, sample_test_logic_Q1) +cat ("sample_test_logic : Me = ") ; cat(q_fig[3]) ; cat (" | seuil = ") ; cat(seuil_logic_Q2) ; cat("\n") ; analyze.perf(res_logic_Q2, "Year", seuil_logic_Q2, sample_test_logic_Q2) +cat ("sample_test_logic : Q3 = ") ; cat(q_fig[4]) ; cat (" | seuil = ") ; cat(seuil_logic_Q3) ; cat("\n") ; analyze.perf(res_logic_Q3, "Year", seuil_logic_Q3, sample_test_logic_Q3) + +######################################################################################## +# Réalisation des prédictions selon un seuil pertinent +# Modèle multimodal sur données d'entrainement + +sample_train_categ_2 <- ifelse(sample_train$Year < q_fig[2], 1, + ifelse(sample_train$Year < q_fig[3], 2, + ifelse(sample_train$Year < q_fig[4], 3, 4))) + +sample_train_categ_predict_2 <- ifelse(res_logic_Q2$fitted.values <= seuil_logic_Q2, + ifelse(res_logic_Q1$fitted.values <= seuil_logic_Q1, 1, 2), + ifelse(res_logic_Q3$fitted.values <= seuil_logic_Q3, 3, 4)) + +cat ("sample_train_categ | seuils pertinents") ; cat("\n") ; analyze.perf.multi(sample_train_categ_2,sample_train_categ_predict_2) + + +# Réalisation des prédictions selon un seuil pertinent +# Modèle multimodal sur données de test + +sample_test_categ_2 <- ifelse(sample_test$Year < q_fig[2], 1, + ifelse(sample_test$Year < q_fig[3], 2, + ifelse(sample_test$Year < q_fig[4], 3, 4))) + +predict_test_Q1 <- predict.glm(res_logic_Q1, sample_test_logic_Q1, type="response") +predict_test_Q2 <- predict.glm(res_logic_Q2, sample_test_logic_Q2, type="response") +predict_test_Q3 <- predict.glm(res_logic_Q3, sample_test_logic_Q3, type="response") + +sample_test_categ_predict_2 <- ifelse(predict_test_Q2 <= seuil_logic_Q2, + ifelse(predict_test_Q1 <= seuil_logic_Q1, 1, 2), + ifelse(predict_test_Q3 <= seuil_logic_Q3, 3, 4)) + +cat ("sample_test_categ | seuils pertinents") ; cat("\n") ; analyze.perf.multi(sample_test_categ_2,sample_test_categ_predict_2) diff --git a/images/Rplot.png b/images/Rplot.png new file mode 100755 index 0000000000000000000000000000000000000000..e34c2511962d5137e42e4ac82ddb511450c74442 Binary files /dev/null and b/images/Rplot.png differ diff --git a/images/Rplot01.png b/images/Rplot01.png new file mode 100755 index 0000000000000000000000000000000000000000..3d8bbcedc187300930d4dbdafe69bf9d1aeacfbd Binary files /dev/null and b/images/Rplot01.png differ diff --git a/images/Rplot02.png b/images/Rplot02.png new file mode 100755 index 0000000000000000000000000000000000000000..89d265eb346b1e94c2a5978210abeb0615804490 Binary files /dev/null and b/images/Rplot02.png differ diff --git a/images/Rplot03.png b/images/Rplot03.png new file mode 100755 index 0000000000000000000000000000000000000000..87272531f1ec3e9963d579a13c8ae7a0a954f68e Binary files /dev/null and b/images/Rplot03.png differ diff --git a/images/Rplot04.png b/images/Rplot04.png new file mode 100755 index 0000000000000000000000000000000000000000..599349825ed654af20b17ed9f33a7afbaf831a62 Binary files /dev/null and b/images/Rplot04.png differ diff --git a/images/Rplot05.png b/images/Rplot05.png new file mode 100755 index 0000000000000000000000000000000000000000..3fffd0355ad66076be116a32f39a80f8814bad22 Binary files /dev/null and b/images/Rplot05.png differ