Skip to content
Snippets Groups Projects
Commit 7e00068c authored by Heloise Chevalier's avatar Heloise Chevalier
Browse files

rapport presque fini

parent b64018c5
No related branches found
No related tags found
No related merge requests found
regresseur <- function(dataset) {
# Chargement de l'environnement
load("env.Rdata")
model.reg <- lm(reg.train.y ~ ., data = reg.train.X)
data.test <- dataset[,reg.mask]
predictions <- predict(model.reg, newdata = data.test)
return(predictions)
}
library(MASS)
library(leaps)
library(corrplot)
library(glmnet)
reg.data <- read.table("tp3_a18_reg_app.txt")
reg.fit <-regsubsets(y ~ .,data=reg.data,method="exhaustive",nvmax=50, really.big = T)
subsets <- summary(reg.fit)
mask <- subsets$which[,2:51]
reg.train.X <- reg.data[,mask]
reg.train.y <- reg.data[,y]
\ No newline at end of file
--- ---
title: "Projet 1 SY19" title: "Projet 1 SY19"
subtitle: "CHEVALIER Héloïse, JORANDON Guillaume" subtitle: "CHEVALIER Héloïse, JORANDON Guillaume"
output: html_notebook output: html_notebook
--- ---
...@@ -17,3 +17,344 @@ Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by press ...@@ -17,3 +17,344 @@ Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by press
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Ctrl+Shift+K* to preview the HTML file). When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Ctrl+Shift+K* to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed. The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
1. Régression
1.1. Analyse exploratoire
Le dataset contient 50 prédicteurs (X1, X2, ..., X50), et une colonne y. Le fichier contient 200 lignes, et toutes les variables sont quantitatives.
On remarque que les prédicteurs sont très peu corrélés entre eux, il ne semble donc pas nécessaire de réaliser une ACP.
```{r}
reg.data <- read.table("tp3_a18_reg_app.txt")
reg.cor <- cor(reg.data)
corrplot(reg.cor)
```
1.2. Recherche du meilleur modèle
Pour rechercher le meilleur modèle, on peut commencer par réaliser une régression linéaire sur l'ensemble du jeu de données, afin de repérer les prédicteurs les plus significatifs.
```{r}
reg <- lm(y ~., data = reg.data)
summary(reg)
```
On compte 13 prédicteurs avec une p-value < 0.001, et 18 prédicteurs avec une p-value < 0.01.
Etant données les nombreuses possibilités de sous-ensembles pertinents pour la régression, on décide d'effecuer une étude exhaustive : on réalise des régressions sur différents sous-ensembles de 2 à 49 prédicteurs, sélectionné avec la méthode pas à pas exhaustive, puis on calcule par cross-validation l'espérance de l'erreur quadratique de chaque modèle.
```{r}
reg.fit <-regsubsets(y ~ .,data=reg.data,method="exhaustive",nvmax=50, really.big = T)
reg.subsets <- summary(reg.fit)
plot(reg.fit,scale="r2", main ="")
title(main = "Choix des meilleurs sous-ensembles par sélection pas à pas exhaustive" )
```
On calcule ainsi l'espérance de l'erreur quadratique de la régression linéaire pour des sous-ensembles de 2 à 48 prédicteurs. On calcule de même l'espérance de l'erreur quadratique pour les régressions ridge et lasso à partir des mêmes sous-ensembles.
La cross-validation se fait en découpant le dataset en 5 parties égales : on calcule l'erreur de test en prennant chacun de ces sous-ensembles comme ensemble de test et les autres comme ensemble d'apprentissage, et on calcule l'erreur quadratique. On répète cette opération 100 fois, puis on calcule la moyenne de toutes les erreurs quadratique. On choisit de faire 100 répétitions pour avoir une moyenne qui ne varie pas trop selon les découpages aléatoires de la cross-validation. En effet, on observe que pour 10 ou 50 répétitions, une valeur d'erreur quadratique peut varier d'une exécution à une autre.
On obtient les courbes suivantes :
```{r}
n.reg = dim(reg.data)[1]
p.reg = dim(reg.data)[2]
reg.models.list <- c("lm", "ridge", "lasso") # Liste des modeles
reg.trials <- 100 # Nombre d'essais de cross-validation par modèle
reg.n_folds <- 5 # Nombre de blocs pour la cross-validation
reg.subset.error <- c()
reg.models.error <- matrix(0,nrow = 48,ncol = 3)
colnames(reg.models.error) <- reg.models.list
rownames(reg.models.error) <- c(2:49)
reg.error.model <- function(model)
{
for (subset in c(2:49))
{
reg.subset.error[subset] <- 0
for(i in c(1:reg.trials)) # On fait plusieurs essais
{
reg.folds <- sample(rep(1:reg.n_folds, length.out=n.reg))
for(k in 1:reg.n_folds) # Sur chaque segment de la cross-validation
{
reg.train <- reg.data[reg.folds!=k,]
reg.train.X <- reg.train[,-p.reg]
reg.train.y <- reg.train$y
reg.test <- reg.data[reg.folds==k,]
reg.test.X <- reg.test[,-p.reg]
reg.test.y <- reg.test$y
reg.subset <- reg.train.X[,reg.subsets$which[subset,2:51]]
reg.subset.test <- reg.test.X[,reg.subsets$which[subset,2:51]]
# Apprentissage du modele
reg.model <- NULL
if(model == "lm")
{
reg.model <- lm(reg.train.y ~., data = reg.subset)
}
else
{
if(model == "ridge")
{
cv.out<-cv.glmnet(as.matrix(reg.subset),reg.train.y,alpha=0)
reg.model<-glmnet(as.matrix(reg.subset),reg.train.y,lambda=cv.out$lambda.min,alpha=0)
}
else if(model == "lasso")
{
cv.out<-cv.glmnet(as.matrix(reg.subset),reg.train.y,alpha=1)
reg.model<-glmnet(as.matrix(reg.subset),reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else
{
reg.model <- NULL
}
}
if(!is.null(reg.model))
{
if(model == "lm")
{
reg.predict <- predict(reg.model, newdata = reg.subset.test)
reg.subset.error[subset] <- reg.subset.error[subset] + mean((reg.test.y - reg.predict)^2)
}
else
{
reg.predict <- predict(reg.model,s=cv.out$lambda.min,newx=as.matrix(reg.subset.test))
reg.subset.error[subset] <- reg.subset.error[subset] + mean((reg.test.y - reg.predict)^2)
}
}
}
}
reg.subset.error[subset] <- reg.subset.error[subset]/(reg.n_folds*reg.trials)
}
reg.subset.error
}
for (model in reg.models.list)
{
reg.models.error[,model] <- as.vector(reg.error.model(model)[2:49])
}
reg.models.error
reg.subsets$which[23,]
```
```{r}
plot(x=c(10:48),y=reg.models.error[10:48,"lm"], col = "red", type = "b", xlab = "nombre de prédicteurs", ylab = "espérance de l'erreur quadratique")
lines(x=c(10:48),y=reg.models.error[10:48,"ridge"], col = "blue", type = "b")
lines(x=c(10:48),y=reg.models.error[10:48,"lasso"], col = "green", type = "b")
```
On observe que l'espérance de l'erreur quadratique la plus faible est obtenue pour 23 prédicteurs, avec la régression lasso, avec une espérance égale à 1271.033.
Il s'agit du sous-ensemble contenant les prédicteurs X1, X2, X3, X10, X14, X16, X17, X18, X19, X24, X28, X31, X32, X34, X35, X37, X38, X39, X40, X41, X43, X46 et X49.
Ces variables ont une p-value comprise entre 2e-16 et 0.12519.
On peut cependant remarquer que cette méthode de sélection d'un modèle de régression n'est pas optimale. En effet, parcourir tous les sous-ensembles en cross-validation est long, et malgré le nombre de répétitions, on observe toujours des variations des valeurs de l'espérance d'une exécution à une autre.
On décide donc de réaliser d'autres sélections de modèles, par la sélection pas à pas avec critères AIC et BIC, ainsi qu'en observant les prédicteurs les plus significatifs.
```{r}
plot(reg.fit,scale="adjr2", main = "")
title(main = "AIC")
plot(reg.fit,scale="bic", main = "")
title(main = "BIC")
```
Le meilleur sous-ensemble selon le critère AIC contient 33 prédicteurs, le sous-ensemble BIC en contient 20.
On peut donc réaliser les régressions linéaire, ridge et lasso sur 4 sous-ensembles : les meilleurs sous-ensembles donnés par les sélections AIC et BIC, et les sous-ensembles des 13 ou 18 prédicteurs les plus significatifs.
On observe sur les courbes présentées précédemment que la régression ridge donne de moins bons résultats que les régressions linéaire et lasso, on décide donc de ne pas effectuer cette régression.
De plus, afin de pouvoir comparer les résultats de tous les modèles pertinents pour les mêmes ensembles d'apprentissage et de test, on recalcule dans la même fonction l'espérance de l'erreur quadratique pour la régression linéaire et la régression lasso des sous-ensembles à 23, 24, 25 et 26 prédicteurs, puisque ce sont ceux pour lesquels les valeurs de l'espérance étaient les plus faibles.
```{r}
reg.models.subsets <- c("lm13", "lasso13", "lm18", "lasso18", "lmAIC", "lassoAIC", "lmBIC", "lassoBIC", "lm23", "lasso23", "lm24", "lasso24", "lm25", "lasso25", "lm26", "lasso26") # Liste des modeles
models.subsets.error <- matrix(0,nrow = 1, ncol = 16)
rownames(models.subsets.error) <- c("espérance")
colnames(models.subsets.error) <- reg.models.subsets
for(model in reg.models.subsets) # Pour chaque modele
{
models.subsets.error[1,model] <- 0
for(i in c(1:reg.trials)) # On fait plusieurs essais
{
reg.folds <- sample(rep(1:reg.n_folds, length.out=n.reg))
for(k in 1:reg.n_folds) # Sur chaque segment de la cross-validation
{
reg.train <- reg.data[reg.folds!=k,]
reg.train.X <- reg.train[,-p.reg]
reg.train.y <- reg.train$y
reg.test <- reg.data[reg.folds==k,]
reg.test.X <- reg.test[,-p.reg]
reg.test.y <- reg.test$y
reg.X.23 <- reg.train.X[,reg.subsets$which[23, 2:51]]
reg.X.24 <- reg.train.X[,reg.subsets$which[24, 2:51]]
reg.X.25 <- reg.train.X[,reg.subsets$which[25, 2:51]]
reg.X.26 <- reg.train.X[,reg.subsets$which[26, 2:51]]
reg.test.23 <- reg.test.X[,reg.subsets$which[23, 2:51]]
reg.test.24 <- reg.test.X[,reg.subsets$which[24, 2:51]]
reg.test.25 <- reg.test.X[,reg.subsets$which[25, 2:51]]
reg.test.26 <- reg.test.X[,reg.subsets$which[26, 2:51]]
# Apprentissage du modele
if(model == "lm13")
{
reg.model <- lm(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 +X41 + X49, data = reg.train)
}
else if(model == "lm18")
{
reg.model <- lm(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 +X41 + X49 + X10 + X24 + X40 + X43 + X46, data = reg.train)
}
else if(model == "lmAIC")
{
reg.model <- lm(y ~ X1 + X2 + X3 + X6 + X9 + X10 + X11 + X14 + X15 + X16 + X17 + X18 + X19 + X21 + X24 + X26 + X27 + X28 + X31 + X32 + X34 + X35 + X37 + X38 + X39 + X40 + X41 + X42 + X43 + X45 + X46 + X47 + X49, data = reg.train)
}
else if(model == "lmBIC")
{
reg.model <- lm(y ~ X1 + X2 + X3 + X10 + X14 + X17 + X19 + X24 + X31 + X32 + X34 + X35 + X37 + X38 + X39 + X40 + X41 + X43 + X46 + X49 , data = reg.train)
}
else if(model == "lm23")
{
reg.model <- lm(reg.train.y ~ . , data = reg.X.23)
}
else if(model == "lm24")
{
reg.model <- lm(reg.train.y ~ . , data = reg.X.24)
}
else if(model == "lm25")
{
reg.model <- lm(reg.train.y ~ . , data = reg.X.25)
}
else if(model == "lm26")
{
reg.model <- lm(reg.train.y ~ . , data = reg.X.26)
}
else if(model == "lasso13")
{
x <- model.matrix(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 +X41 + X49, reg.data)
xapp <- x[reg.folds!=k,]
xtst <- x[reg.folds == k,]
cv.out<-cv.glmnet(xapp,reg.train.y,alpha=1)
reg.model<-glmnet(xapp,reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else if(model == "lasso18")
{
x <- model.matrix(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 +X41 + X49 + X10 + X24 + X40 + X43 + X46, reg.data)
xapp <- x[reg.folds!=k,]
xtst <- x[reg.folds == k,]
cv.out<-cv.glmnet(xapp,reg.train.y,alpha=1)
reg.model<-glmnet(xapp,reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else if(model == "lassoAIC")
{
x <- model.matrix(y ~ X1 + X2 + X3 + X6 + X9 + X10 + X11 + X14 + X15 + X16 + X17 + X18 + X19 + X21 + X24 + X26 + X27 + X28 + X31 + X32 + X34 + X35 + X37 + X38 + X39 + X40 + X41 + X42 + X43 + X45 + X46 + X47 + X49, reg.data)
xapp <- x[reg.folds!=k,]
xtst <- x[reg.folds == k,]
cv.out<-cv.glmnet(xapp,reg.train.y,alpha=1)
reg.model<-glmnet(xapp,reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else if(model == "lassoBIC")
{
x <- model.matrix(y ~ X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 + X1 + X49, reg.data)
xapp <- x[reg.folds!=k,]
xtst <- x[reg.folds == k,]
cv.out<-cv.glmnet(xapp,reg.train.y,alpha=1)
reg.model<-glmnet(xapp,reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else if(model == "lasso23")
{
cv.out<-cv.glmnet(as.matrix(reg.X.23),reg.train.y,alpha=1)
reg.model<-glmnet(as.matrix(reg.X.23),reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else if(model == "lasso24")
{
cv.out<-cv.glmnet(as.matrix(reg.X.24),reg.train.y,alpha=1)
reg.model<-glmnet(as.matrix(reg.X.24),reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else if(model == "lasso25")
{
cv.out<-cv.glmnet(as.matrix(reg.X.25),reg.train.y,alpha=1)
reg.model<-glmnet(as.matrix(reg.X.25),reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else if(model == "lasso26")
{
cv.out<-cv.glmnet(as.matrix(reg.X.26),reg.train.y,alpha=1)
reg.model<-glmnet(as.matrix(reg.X.26),reg.train.y,lambda=cv.out$lambda.min,alpha=1)
}
else
{
reg.model <- NULL
}
if(!is.null(reg.model))
{
if(model == "lm13" || model == "lm18" || model == "lmAIC" || model == "lmBIC" || model == "lm23" || model == "lm24" || model == "lm25" || model == "lm26")
{
reg.predict <- predict(reg.model, newdata = reg.test)
models.subsets.error[1,model] <- models.subsets.error[1,model] + mean((reg.test.y - reg.predict)^2)
}
else if (model == "lasso13" || model == "lasso18" || model == "lassoAIC" || model == "lassoBIC")
{
reg.predict <- predict(reg.model,s=cv.out$lambda.min,newx=xtst)
models.subsets.error[1,model] <- models.subsets.error[1,model] + mean((reg.test.y - reg.predict)^2)
}
else if (model == "lasso23")
{
reg.predict <- predict(reg.model,s=cv.out$lambda.min,newx=as.matrix(reg.test.23))
models.subsets.error[1,model] <- models.subsets.error[1,model] + mean((reg.test.y - reg.predict)^2)
}
else if (model == "lasso24")
{
reg.predict <- predict(reg.model,s=cv.out$lambda.min,newx=as.matrix(reg.test.24))
models.subsets.error[1,model] <- models.subsets.error[1,model] + mean((reg.test.y - reg.predict)^2)
}
else if (model == "lasso25")
{
reg.predict <- predict(reg.model,s=cv.out$lambda.min,newx=as.matrix(reg.test.25))
models.subsets.error[1,model] <- models.subsets.error[1,model] + mean((reg.test.y - reg.predict)^2)
}
else if (model == "lasso26")
{
reg.predict <- predict(reg.model,s=cv.out$lambda.min,newx=as.matrix(reg.test.26))
models.subsets.error[1,model] <- models.subsets.error[1,model] + mean((reg.test.y - reg.predict)^2)
}
}
}
}
models.subsets.error[1,model] <- models.subsets.error[1,model]/(reg.n_folds*reg.trials)
}
print(models.subsets.error)
```
```{r}
plot(x = c(1:16), y = models.subsets.error, xlab = "modèle", ylab = "espérance de l'erreur quadratique")
```
On remarque que la plus faible valeur de l'espérance est obtenue cette fois pour le modèle de régression linéaire avec un sous-ensemble de 23 prédicteurs. On observe cependant que le modèle de régression lasso avec 23 prédicteurs obtient également une espérance parmis les plus faibles valeurs, ainsi que les régressions linéaire et lasso à 24 prédicteurs.
On choisit donc le modèle suivant : régression linéaire avec les 23 prédicteurs : X1, X2, X3, X10, X14, X16, X17, X18, X19, X24, X28, X31, X32, X34, X35, X37, X38, X39, X40, X41, X43, X46 et X49.
1.3. Analyse du modèle sélectionné
Q-Q plot entre -2 et 2
Distance de Cook pour les points aberrants
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment