Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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
#sur toutes les donnes du coup
reg <- lm(y ~ ., data = reg.data)
abline(0,1)
#prdicteurs les plus significatifs:
# X1, X2, X3, -X10, X14, X19, -X24, X32, X34, X35, X37, X38, X39, -X40, X41, -X43, -X46, X49
#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))
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#prdicteurs 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)
#esprance 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)
#esprance 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 donnes)
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 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 + X24 + X46 + X49, data = reg.data)
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
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
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
#meilleur resultat so far (pas de beaucoup par rapport au modle 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