bootstrap sur une régression multiple

Version imprimable

Bonjour à tous,

Je souhaite réaliser un bootstrap sur une régression linéaire multiple. En effet, je souhaite tester plusieurs fois mes modèles avec à chaque répétition un nouveau jeu d'apprentissage et de test pris au hasard d'un jeu de donnée initial (ici d). les modèles sont crées à partir de la fonction "regsubsets" du package "leaps".

voici mon code:
Code:

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 39 40 41 42 43 44 45 46 47 n<-100 d<-data.frame (y=runif(n), x1=(runif(n))*0.5, x2=(runif(n))*0.25, x3=(runif (n))*3.2, x4=(runif(n))*4.35, x5=(runif(n))*2.2) train<-sample(seq(100), 70, replace=FALSE) regfit.bestfwd=regsubsets(y~.,data=d[train,], nvmax =5, method="forward") summar.regfitfwd<-summary(regfit.bestfwd) #le calcul de la statistique des modèles sur un jeu d'apprentissage app.mat=model.matrix(y~.,data=d[train,]) {val.errors=rep(NA,5) stor.r2=rep(NA,5) stor.rmse=rep(NA,5) stor.pe=rep(NA,5) for(i in 1:5){ coefi=coef(regfit.bestfwd,id=i) pred=app.mat [,names(coefi)]%*% coefi val.errors[i]=mean((d$y[train]-pred)^2) stor.r2[i]=1-(sum((d$y[train]-pred)^2)/sum((d$y[train]-mean(d$y[train]))^2)) stor.rmse[i]=(sum((pred-d$y[train])^2/20))^0.5 #rmse: root mean squared error stor.pe[i]=stor.rmse[i]/mean(d$y[train]) #pe: pourcentage d'erreur }} stor.r2 which.max(summar.regfitfwd$rsq) coef(regfit.bestfwd, 5) #le calcul de la statistique des modèles sur un jeu test app.mat=model.matrix(y~.,data=d[-train,]) {val.errors=rep(NA,5) stor.r2=rep(NA,5) stor.rmse=rep(NA,5) stor.pe=rep(NA,5) for(i in 1:5){ coefi=coef(regfit.bestfwd,id=i) pred=app.mat [,names(coefi)]%*% coefi val.errors[i]=mean((d$y[-train]-pred)^2) stor.r2[i]=1-(sum((d$y[-train]-pred)^2)/sum((d$y[-train]-mean(d$y[-train]))^2)) stor.rmse[i]=(sum((pred-d$y[-train])^2/20))^0.5 stor.pe[i]=stor.rmse[i]/mean(d$y[-train]) }} stor.r2 which.max(summar.regfitfwd$rsq) coef(regfit.bestfwd, 5)
Au fait je souhaite répéter autant de fois que possible le code ci-dessus avec à chaque répétition un nouveau jeu de donnée d'apprentissage et de test pris au hasard à partir du jeu (d). Cela me permettra d'obtenir la statistique descriptive (moyenne, écart-type) du R², rmse et pe du modèle à une variable descriptive puis à deux variable descriptives puis à trois...etc jusqu'à cinq variables.

J'espère que c'est clair

Merci.

Bonjour,

oui c'est clair, mais tu parles de bootstrap mais ça n'en est pas un. Un bootstrap tu tires au hasard n individus parmi n et avec remise.
Ici si je comprends bien ton code, tu sépares ton jeu de données en 70% pour la calibration et 30% pour le test.
Pour le reste il te suffit de mettre tout ton code dans une boucle sur j et de créer des objets de stockage à plusieurs dimension comme des matrices ou des array.

Code:

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
39
40
41
42
 
n<-100
d<-data.frame (y=runif(n), x1=(runif(n))*0.5, x2=(runif(n))*0.25, x3=(runif (n))*3.2,
               x4=(runif(n))*4.35, x5=(runif(n))*2.2)
 
N <- 999 # le nombre de répétition
 
val.errorst=matrix(numeric(),N, 5)
stor.r2t=matrix(numeric(),N, 5)
stor.rmset=matrix(numeric(),N, 5)
stor.pet=matrix(numeric(),N, 5)
 
val.errors=matrix(numeric(),N, 5)
stor.r2=matrix(numeric(),N, 5)
stor.rmse=matrix(numeric(),N, 5)
stor.pe=matrix(numeric(),N, 5)
 
for (j in 1:N) {
  train<-sample(seq(100), 70, replace=FALSE)
  regfit.bestfwd=regsubsets(y~.,data=d[train,],
                            nvmax =5, method="forward")
  summar.regfitfwd<-summary(regfit.bestfwd)
 
  #le calcul de la statistique des modèles sur un jeu d'apprentissage 
  app.matt = model.matrix(y~.,data=d[train,])
  #le calcul de la statistique des modèles sur un jeu test 
  app.mat=model.matrix(y~.,data=d[-train,])
  for(i in 1:5){
    coefi=coef(regfit.bestfwd,id=i)
    pred=app.matt[,names(coefi)]%*% coefi
    val.errorst[j, i]=mean((d$y[train]-pred)^2)
    stor.r2t[j, i]=1-(sum((d$y[train]-pred)^2)/sum((d$y[train]-mean(d$y[train]))^2))
    stor.rmset[j, i]=(sum((pred-d$y[train])^2/20))^0.5 #rmse: root mean squared error
    stor.pet[j, i]=stor.rmset[j, i]/mean(d$y[train]) #pe: pourcentage d'erreur
 
    pred=app.mat [,names(coefi)]%*% coefi
    val.errors[j, i]=mean((d$y[-train]-pred)^2)
    stor.r2[j, i]=1-(sum((d$y[-train]-pred)^2)/sum((d$y[-train]-mean(d$y[-train]))^2))
    stor.rmse[j, i]=(sum((pred-d$y[-train])^2/20))^0.5
    stor.pe[j, i]=stor.rmse[j, i]/mean(d$y[-train])
  }
}

Après il faut que tu te méfies. A chaque tirage aléatoire tu vas peut-être avoir des sélections de variables différentes pour les modèles avec 1 à p-1 prédicteurs. Si tel est le cas ça peut revenir a comparer des choux et des carottes, tout dépend de tes objectifs.

04/10/2018, 11h22
Lou.lou

Merci infiniment Tototod, je vois que tu as consacré du temps à réécrire le code, c'est très simpat de ta part. Vraiment merci.

Citation:

Un bootstrap tu tires au hasard n individus parmi n et avec remise

Merci pour le rappel, du coup pour que ça devienne un bootstrap il suffit de mettre dans le code sampl(replace=TRUE).

Ton code marche à merveille, cependant je cale sur le R², apparemment il y a une erreur à ce niveau qui m'échappe, le R² est toujours négatif et<-1. Est-ce ma formule qui est fausse?

Code:

stor.r2[j, i]=1-(sum((d$y[-train]-pred)^2)/sum((d$y[-train]-mean(d$y[-train]))^2))

Citation:

A chaque tirage aléatoire tu vas peut-être avoir des sélections de variables différentes pour les modèles avec 1 à p-1 prédicteurs

Sur la base de ton commentaire ci-dessus, voici ce que je suggère:
1-Sélectionner un pool de modèles sur la base de la moyenne du R² le plus élevé. Exemple les modèles à 3 variables.
2-Dans ce pool de modèles, prendre les 3 variables qui sont le plus souvent sélectionnées ensemble.

Si mon raisonnement est valide, est-ce que la fonction ci-dessous me sélectionnera ces 3 variables :

Code:

coef(regfit.bestfwd,3)

Si non il faudra écrire un code pour le faire, et là je ne vois pas trop comment m'y prendre :ptdr:.
04/10/2018, 16h09
tototode

Concernant la stratégie c'est à toi de voir ce que tu veux faire.
Tu peux aussi partir du résultat sur toutes les données et voir ce que ça donne en bootstrap en figeant les choses dès le début.
Tu devrais aussi regarder plus en détail la méthode boostrap et ce qu'elle peut t'apporter.
Tu peux aussi t'orienter vers de la validation croisée.

Tout dépend tes objectifs.
Cdlt
04/10/2018, 16h53
Lou.lou

Pour ce qui est du R² j'ai utilisé la fonction "(cor(x,y))^2" et ça marche.

attention il y a juste une erreur dans le calcul du rmse où je divise par 20 alors qu'il faut diviser par nrow(d[train]) pour le rmse sur la base de données d'apprentissage et par nrow(d[-train]) pour le rmse sur la base de données tests.

Merci pour les conseils Tototode et à bientôt.