From 5e18fb0beb5a5e00306a9f990f14aebed69aab60 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?H=C3=A9lo=C3=AFse=20Chevalier?=
 <heloise.chevalier@etu.utc.fr>
Date: Tue, 16 Oct 2018 12:13:35 +0200
Subject: [PATCH] ridge & lasso regression

---
 "r\303\251gression.r" | 199 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 199 insertions(+)
 create mode 100644 "r\303\251gression.r"

diff --git "a/r\303\251gression.r" "b/r\303\251gression.r"
new file mode 100644
index 0000000..9d5e8a9
--- /dev/null
+++ "b/r\303\251gression.r"
@@ -0,0 +1,199 @@
+library(MASS)
+library(leaps)
+library(corrplot)
+library(glmnet)
+
+#chargement des données
+reg.data <- read.table("tp3_a18_reg_app.txt")
+
+x.reg <- reg.data[,1:50]
+y.reg <- reg.data[,"y"]
+
+n = dim(reg.data)[1]
+
+reg.mask = sort(sample(x = 1:n, size = trunc(n*2/3)))
+
+reg.appr = reg.data[reg.mask,]
+reg.test = reg.data[-reg.mask,]
+x.appr = reg.appr[,1:50]
+y.appr = reg.appr[,"y"]
+
+x.test = reg.test[,1:50]
+y.test = reg.test[,"y"]
+
+#stocker modèle dans le fichier RData
+#pas stocker tous nos essais dans .Rdata (duh)
+#fonctions : doivent retourner des prédictions à partir du jeu de données en argument
+#d'après Sylvain Rousseau : pour la régression on a pas à calculer d'erreurs de test
+#                           on apprend sur toutes les données
+# mais d'après le sujet il faut calculer l'espérance de l'erreur quadratique donc???
+
+#il n'y aura pas la colonne y dans le data de test des profs
+
+#infos sur données
+summary(reg.data)
+cor <- cor(reg.data)
+corrplot(cor) #prédicteurs pas corrélés entre eux (à part X19 corrélé avec y maybe???)
+
+
+#première régression test pour voir prédicteurs les plus significatifs
+reg <- lm(y ~ ., data = reg.appr)
+summary(reg)
+confint(reg)
+
+plot(y.appr,reg$fitted.values)
+abline(0,1)
+
+#prédicteurs les plus significatifs:
+# X1, X2, X3, -X10, X14, X19, -X24, X32, X34, X35, X37, X38, X39, -X40, X41, -X43
+
+#espérance de l'erreur quadratique :
+pred <- predict(reg,newdata = reg.test)
+plot(y.test,pred)
+abline(0,1)
+
+mean((y.test - pred)^2)
+#2052.946
+
+#infos sur les résidus
+rres = reg$residuals
+rstd = rstandard(reg)
+rstu = rstudent(reg)
+
+plot(y.appr,rstd)
+plot(y.appr,rstu)
+
+shapiro.test(rres)
+
+## Q-Q plots
+qqnorm(rres, asp = 1)
+qqline(rres, dist = qnorm)
+
+qqnorm(rstd, asp = 1)
+qqline(rstd, dist = qnorm)
+
+qqnorm(rstu, asp = 1)
+qqline(rstu, dist = qnorm)
+
+#influence globale
+plot(reg, which = 4, cook.levels = c(0, 0.1))
+plot(reg, which = 5, cook.levels = c(0, 0.1))
+#128 aberrant
+
+
+#stepwise selection
+reg.fit.exh <-regsubsets(y ~ .,data=reg.data,method="exhaustive",nvmax=15, really.big = T)
+x11()
+plot(reg.fit.exh,scale="r2", main ="")
+title(main = "best subset" )
+
+
+#forward
+reg.fit.f<-regsubsets(y ~ .,data=reg.data,method="forward",nvmax=15, really.big = T)
+plot(reg.fit.f,scale="r2", main = "")
+title(main ="forward")
+
+#backward
+reg.fit.b<-regsubsets(y ~ .,data=reg.data,method="backward",nvmax=15, really.big = T)
+plot(reg.fit.b,scale="r2", main = "")
+title(main = "backward")
+
+#AIC et BIC
+reg.fit <- regsubsets(y ~., data = reg.data, method = "exhaustive", nvmax = 15, really.big = T)
+x11()
+plot(reg.fit,scale="adjr2", main = "") 
+title(main = "AIC")
+x11()
+plot(reg.fit,scale="bic", main = "")
+title(main = "BIC")
+
+#pas de différence entre tous ces résultats, on part sur un premier subset avec tous les prédicteurs au dessus de 0.9
+
+reg.model <- lm(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 + X41, data = reg.appr)
+summary(reg.model) #X24, X46 et X49 pas très significatif ici, ils dégagent
+confint(reg.model)
+
+plot(y.appr,reg.model$fitted.values)
+abline(0,1)
+
+#espérance de l'erreur quadratique :
+pred.model <- predict(reg.model,newdata = reg.test)
+plot(y.test,pred.model)
+abline(0,1)
+
+mean((y.test - pred.model)^2)
+#1755.945 c mieux
+
+#infos sur les résidus
+rres.model = reg.model$residuals
+rstd.model = rstandard(reg.model)
+rstu.model = rstudent(reg.model)
+
+plot(y.appr,rstd.model)
+plot(y.appr,rstu.model)
+
+shapiro.test(rres.model)
+
+## 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)
+
+#influence globale
+plot(reg.model, which = 4, cook.levels = c(0, 0.1))
+plot(reg.model, which = 5, cook.levels = c(0, 0.1))
+#aucun point aberrant on est contents (quelques point un peu limites : 114, 118, 140)
+
+#transfo non linéaires ?
+# les plots des résidus en fonction de chaque prédicteur ne révèlent aucune tendance qui suggèrerait une transfo non linéaire
+plot(reg.appr$X1,rstd.model)
+plot(reg.appr$X2,rstd.model)
+plot(reg.appr$X3,rstd.model)
+plot(reg.appr$X14,rstd.model)
+plot(reg.appr$X17,rstd.model)
+plot(reg.appr$X19,rstd.model) #quadratique?
+plot(reg.appr$X32,rstd.model) 
+plot(reg.appr$X34,rstd.model) #quadratique?
+plot(reg.appr$X35,rstd.model)
+plot(reg.appr$X37,rstd.model)
+plot(reg.appr$X38,rstd.model)
+plot(reg.appr$X39,rstd.model)
+plot(reg.data$X41,rstd.model)
+
+#à voir
+
+##ridge regression
+x <- model.matrix(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 + X41, reg.data)
+xapp <- x[reg.mask,]
+xtst <- x[-reg.mask,]
+
+cv.out<-cv.glmnet(xapp,y.appr,alpha=0)
+plot(cv.out)
+
+fit<-glmnet(xapp,y.appr,lambda=cv.out$lambda.min,alpha=0)
+
+#esperance de l'erreur quadratique :
+pred.ridge <- predict(fit,s=cv.out$lambda.min,newx = xtst)
+plot(y.test,pred.ridge)
+abline(0,1)
+
+mean((y.test - pred.ridge)^2)
+#1857.962 on a vu mieux
+
+
+#lasso regression
+cv.out.lasso <- cv.glmnet(xapp,y.appr,alpha=1)
+plot(cv.out.lasso)
+fit.lasso <- glmnet(xapp,y.appr,lambda=cv.out.lasso$lambda.min,alpha=1)
+pred.lasso <- predict(fit.lasso,s=cv.out.lasso$lambda.min,newx=xtst)
+plot(y.test,pred.lasso)
+abline(0,1)
+
+mean((y.test - pred.lasso)^2)
+#1754.978 
+#meilleur resultat so far (pas de beaucoup par rapport au modèle avec subset)
\ No newline at end of file
-- 
GitLab