rapport.Rmd

title: "Projet 1 SY19"
subtitle: CHEVALIER Heloise, JORANDON Guillaume
output: html_notebook

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.

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.

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.

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 :

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])
}

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.

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])
}

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.

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.

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)
}

#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é.

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 :

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")

shapiro.test(rres.model)

Le test de Shapiro appliqué aux résidus donne W = 0.99286, ce qui confirme l'hypothèse de normalité.

## 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é :

#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.