---
title: "Projet 1 SY19"
subtitle: CHEVALIER Heloise, JORANDON Guillaume
output: html_notebook
---
```{r}
library(MASS)
library(leaps)
library(corrplot)
library(glmnet)
```
#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 de ces 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 10 fois dans un premier temps, puis on calcule la moyenne de toutes les erreurs quadratique.
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.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, reg.n_folds, reg.trials, which.subsets)
{
for (subset in which.subsets)
{
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, 5,10, c(2:49))[2:49])
}
```
```{r}
plot(reg.models.error[,"lm"], col = "red", type = "l", xlab = "nombre de pr�dicteurs", ylab = "esp�rance de l'erreur quadratique")
lines(reg.models.error[,"ridge"], col = "blue", type = "l")
lines(reg.models.error[,"lasso"], col = "green", type = "l")
legend(30,7000,legend = c("lin�aire", "ridge", "lasso"), col = c("red","blue","green"), lty = 1:1)
```
On observe que l'esp�rance de l'erreur quadratique en fonction du nombre de pr�dicteurs diminue jusqu'� atteindre un minimum pour un nombre de pr�dicteurs compris dans l'intervalle [15;30].
On r�p�te donc l'op�ration, mais en se limitant aux sous-ensembles de 20 � 26 pr�dicteurs.
On peut alors augmenter le nombre d'essais � 100 pour obtenir des moyennes plus stables : en effet, on observe que pour 10 r�p�titions, l'esp�rance d'un mod�le varie d'une ex�cution � une autre.
```{r}
reg.subset.error <- c()
reg.models.error.zoom <- matrix(0,nrow = 16,ncol = 3)
colnames(reg.models.error.zoom) <- reg.models.list
rownames(reg.models.error.zoom) <- c(15:30)
for (model in reg.models.list)
{
reg.models.error.zoom[,model] <- as.vector(reg.error.model(model, 5,100, c(15:30))[15:30])
}
```
```{r}
plot(x=c(15:30),y=reg.models.error.zoom[,"lm"], col = "red", type = "b", xlab = "nombre de pr�dicteurs", ylab = "esp�rance de l'erreur quadratique")
lines(x=c(15:30),y=reg.models.error.zoom[,"ridge"], col = "blue", type = "b")
lines(x=c(15:30),y=reg.models.error.zoom[,"lasso"], col = "green", type = "b")
legend(23,1425,legend = c("lin�aire", "ridge", "lasso"), col = c("red","blue","green"), lty = 1:1)
```
On remarque �galement que la r�gression ridge obtient globalement de moins bons r�sultats, quel que soit le sous-ensemble de pr�dicteurs.
L'esp�rance de l'erreur quadratique la plus faible est obtenue ici pour 24 pr�dicteurs, avec la r�gression lasso.
Il s'agit du sous-ensemble contenant les pr�dicteurs X1, X2, X3, X9, X10, X14, X16, X17, X18, X19, X24, X27, X28, 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 selon le crit�re 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.
La r�gression ridge donnant de moins bons r�sultats, on d�cide de n'effectuer que les r�gressions lin�aire et lasso.
De plus, afin de pouvoir comparer les r�sultats de tous les mod�les pertinents pour les m�mes ensembles d'apprentissage et de validation, 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.trials <- 100
reg.n_folds <- 5
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 = 16, ncol = 1)
rownames(models.subsets.error) <- reg.models.subsets
colnames(models.subsets.error) <- c("esp�rance")
#print(models.subsets.error["lm13", 1])
for(model in reg.models.subsets) # Pour chaque modele
{
#models.subsets.error[model,1] <- 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[model,1] <- models.subsets.error[model,1] + 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[model,1] <- models.subsets.error[model,1] + 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[model,1] <- models.subsets.error[model,1] + 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[model,1] <- models.subsets.error[model,1] + 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[model,1] <- models.subsets.error[model,1] + 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[model,1] <- models.subsets.error[model,1] + mean((reg.test.y - reg.predict)^2)
}
}
}
}
models.subsets.error[model,1] <- models.subsets.error[model,1]/(reg.n_folds*reg.trials)
}
```
```{r}
#plot(x = c(1:16), y = models.subsets.error, xlab = "mod�le", ylab = "esp�rance de l'erreur quadratique")
```
esp�rance
lm13 1539.681
lasso13 1541.783
lm18 1348.237
lasso18 1345.840
lmAIC 1347.019
lassoAIC 1367.446
lmBIC 1298.314
lassoBIC 1620.786
lm23 1274.732
lasso23 1292.388
lm24 1273.862
lasso24 1266.733
lm25 1276.925
lasso25 1277.970
lm26 1274.081
lasso26 1283.311
On remarque que la plus faible valeur de l'esp�rance est obtenue cette fois encore pour le mod�le de r�gression lasso avec un sous-ensemble de 24 pr�dicteurs.
On peut noter cependant que l'�cart entre les plus faibles valeurs de l'esp�rance est toujours de + ou - 10 (valeurs comprises entre 1266 et 1277), les valeurs minimum �tant obtenues pour les sous-ensembles � 23, 24, 25 et 26 pr�dicteurs. De plus, �tant donn� que les r�sultats de ces mod�les sont proches, et que les valeurs de l'esp�rance varient d'une ex�cution � l'autre, le mod�le obtenant le meilleur r�sultat varie �galement d'une ex�cution � l'autre.
Pour d�terminer le meilleur mod�le parmis ces 8 possibilit�s (r�gressions lin�aire et lasso pour les 4 sous-ensembles obtenant les meilleurs r�sultats), on r�alise une derni�re cross-validation en ne conservant que ces mod�les, mais avec un nombre d'essai plus �lev�.
```{r}
reg.trials <- 250
reg.n_folds <- 5
reg.models.subsets.best <- c("lm23", "lasso23", "lm24", "lasso24", "lm25", "lasso25", "lm26", "lasso26") # Liste des modeles
models.subsets.error.best <- matrix(0,nrow = 8, ncol = 1)
rownames(models.subsets.error.best) <- reg.models.subsets.best
colnames(models.subsets.error.best) <- c("esp�rance")
for(model in reg.models.subsets.best) # Pour chaque modele
{
#models.subsets.error.best[model,1] <- 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 == "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 == "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 == "lm23" || model == "lm24" || model == "lm25" || model == "lm26")
{
reg.predict <- predict(reg.model, newdata = reg.test)
models.subsets.error.best[model,1] <- models.subsets.error.best[model,1] + 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.best[model,1] <- models.subsets.error.best[model,1] + 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.best[model,1] <- models.subsets.error.best[model,1] + 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.best[model,1] <- models.subsets.error.best[model,1] + 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.best[model,1] <- models.subsets.error.best[model,1] + mean((reg.test.y - reg.predict)^2)
}
}
}
}
models.subsets.error.best[model,1] <- models.subsets.error.best[model,1]/(reg.n_folds*reg.trials)
}
print(models.subsets.error.best)
```
esp�rance
lm23 1275.108
lasso23 1273.076
lm24 1264.330
lasso24 1269.377
lm25 1274.778
lasso25 1271.424
lm26 1280.144
lasso26 1280.023
Le meilleur mod�le est ici la r�gression lin�aire � 24 pr�dicteurs, bien que l'esp�rance de la r�gression lasso avec les m�mes pr�dicteurs soit tr�s proche.
On choisit donc le mod�le suivant : r�gression lin�aire avec les 24 pr�dicteurs : X1, X2, X3, X9, X10, X14, X16, X17, X18, X19, X24, X27, X28, X32, X34, X35, X37, X38, X39, X40, X41, X43, X46 et X49.
##1.3. Analyse du mod�le s�lectionn�
###Analyse des r�sidus :
```{r}
model.fit <- lm(y ~ X1 + X2 + X3 + X9 + X10 + X14 + X16 + X17 + X18 + X19 + X24 + X27 + X28 + X32 + X34 + X35 + X37 + X38 + X39 + X40 + X41 + X43 + X46 + X49, data = reg.data)
rres.model = model.fit$residuals
rstd.model = rstandard(model.fit)
rstu.model = rstudent(model.fit)
plot(reg.data$y,rstd.model, main ="")
title("R�sidus standaris�s en fonction de y")
plot(reg.data$y,rstu.model, main = "")
title("R�sidus studentis�s en fonction de y")
```
```{r}
shapiro.test(rres.model)
```
Le test de Shapiro appliqu� aux r�sidus donne W = 0.99286, ce qui confirme l'hypoth�se de normalit�.
```{r}
## Q-Q plots
qqnorm(rres.model, asp = 1)
qqline(rres.model, dist = qnorm)
qqnorm(rstd.model, asp = 1)
qqline(rstd.model, dist = qnorm)
qqnorm(rstu.model, asp = 1)
qqline(rstu.model, dist = qnorm)
```
De m�me, les Q-Q plots des r�sidus standardis�s et studentis�s montrent que la plupart des valeurs sont comprises entre -2 et 2.
On pourrait �galement envisager une r�gression non lin�aire, mais l'analyse des r�sidus en fonction de chacun des pr�dicteurs du mod�le ne montre aucune variation du bruit en fonction d'un pr�dicteur, ce qui ne nous incite pas � complexifier le mod�le.
###Analyse de la stabilit� :
```{r}
#influence globale
plot(model.fit, which = 4, cook.levels = c(0, 0.1))
plot(model.fit, which = 5, cook.levels = c(0, 0.1))
```
L'utilisation de la distance de Cook ne montre aucune valeur aberrante, bien que le point 118 se d�marque des autres par une valeur de distance plus �lev�e.
La distance de Cook du point 118 �tant < 0.1, il ne semble cependant pas n�cessaire de supprimer cette ligne.