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 
#sur toutes les donn�es du coup
reg <- lm(y ~ ., data = reg.data)
summary(reg)
confint(reg)

plot(reg.data$y,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, -X46, X49

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

#pr�dicteurs les plus significatifs:
# X1, X2, X3, -X10, X14, X19, -X24, X32, X34, X35, X37, X38, X39, -X40, X41, -X43, -X46, X49

model.13pred <- lm(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 +X41 + X49, data = reg.appr)
summary(model.13pred) 
confint(model.13pred)

plot(y.appr,model.13pred$fitted.values)
abline(0,1)

#esp�rance de l'erreur quadratique :
pred.model.13 <- predict(model.13pred,newdata = reg.test)
plot(y.test,pred.model.13)
abline(0,1)

mean((y.test - pred.model.13)^2)
#1792.472 euh ouais bof

model.18pred <- lm(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 +X41 + X49 + X10 + X24 + X40 + X43 + X46, data = reg.appr)
summary(model.18pred) 
confint(model.18pred)

plot(y.appr,model.18pred$fitted.values)
abline(0,1)

#esp�rance de l'erreur quadratique :
pred.model.18 <- predict(model.18pred,newdata = reg.test)
plot(y.test,pred.model.18)
abline(0,1)

mean((y.test - pred.model.18)^2)
#1730.831 


#stepwise selection (sur toutes les donn�es)
reg.fit.exh <-regsubsets(y ~ .,data=reg.data,method="exhaustive",nvmax=22, 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=22, really.big = T)
X11()
plot(reg.fit.b,scale="r2", main = "")
title(main = "backward")

#AIC et BIC
reg.fit <- regsubsets(y ~., data = reg.data, method = "exhaustive", nvmax = 18, 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 + X24 + X46 + X49, data = reg.data)
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)




##ridge regression
x18 <- model.matrix(y ~ X1 + X2 + X3 + X14 + X19 + X32 + X34 + X35 + X37 + X38 + X39 +X41 + X49 + X10 + X24 + X40 + X43 + X46, reg.data)
xapp18 <- x18[reg.mask,]
xtst18 <- x18[-reg.mask,]

cv.out18<-cv.glmnet(xapp18,y.appr,alpha=0)
plot(cv.out18)

fit18<-glmnet(xapp18,y.appr,lambda=cv.out18$lambda.min,alpha=0)

#esperance de l'erreur quadratique :
pred.ridge18 <- predict(fit18,s=cv.out18$lambda.min,newx = xtst18)
plot(y.test,pred.ridge18)
abline(0,1)

mean((y.test - pred.ridge18)^2)
#1699.158


#lasso regression
cv.out.lasso18 <- cv.glmnet(xapp18,y.appr,alpha=1)
plot(cv.out.lasso18)
fit.lasso18 <- glmnet(xapp18,y.appr,lambda=cv.out.lasso18$lambda.min,alpha=1)
pred.lasso18 <- predict(fit.lasso18,s=cv.out.lasso18$lambda.min,newx=xtst18)
plot(y.test,pred.lasso18)
abline(0,1)

mean((y.test - pred.lasso18)^2)
#1707.968