ridge & lasso regression

régression.r

+library(MASS)
+library(leaps)
+library(corrplot)
+library(glmnet)
+
+#chargement des donnes
+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 modle dans le fichier RData
+#pas stocker tous nos essais dans .Rdata (duh)
+#fonctions : doivent retourner des prdictions  partir du jeu de donnes en argument
+#d'aprs Sylvain Rousseau : pour la rgression on a pas  calculer d'erreurs de test
+#                           on apprend sur toutes les donnes
+# mais d'aprs le sujet il faut calculer l'esprance de l'erreur quadratique donc???
+
+#il n'y aura pas la colonne y dans le data de test des profs
+
+#infos sur donnes
+summary(reg.data)
+cor <- cor(reg.data)
+corrplot(cor) #prdicteurs pas corrls entre eux ( part X19 corrl avec y maybe???)
+
+
+#premire rgression test pour voir prdicteurs les plus significatifs
+reg <- lm(y ~ ., data = reg.appr)
+summary(reg)
+confint(reg)
+
+plot(y.appr,reg$fitted.values)
+abline(0,1)
+
+#prdicteurs les plus significatifs:
+# X1, X2, X3, -X10, X14, X19, -X24, X32, X34, X35, X37, X38, X39, -X40, X41, -X43
+
+#esprance 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 rsidus
+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 diffrence entre tous ces rsultats, on part sur un premier subset avec tous les prdicteurs 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 trs significatif ici, ils dgagent
+confint(reg.model)
+
+plot(y.appr,reg.model$fitted.values)
+abline(0,1)
+
+#esprance 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 rsidus
+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 linaires ?
+# les plots des rsidus en fonction de chaque prdicteur ne rvlent aucune tendance qui suggrerait une transfo non linaire
+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 modle avec subset)
\ No newline at end of file