Contexte
On s’intéresse au lien éventuel entre le poids d’un homme et divers caractéristiques physiques. Pour 22 hommes en bonne santé âgés de 16 à 30 ans, on dispose :
- du poids en kg (variable \(Y\) ),
- de la circonférence maximale de l’avant-bras en cm (variable \(X_1\)),
- de la circonférence maximale du biceps en cm (variable \(X_2\)),
- de la distance autour de la poitrine directement sous les aisselles en cm (variable \(X_3\)),
- de la distance autour du cou, à peu près à mi-hauteur, en cm (variable \(X_4\)),
- de la distance autour des épaules, mesurées autour de la pointe des omoplates, en cm (variable \(X_5\)),
- de la distance autour de la taille au niveau de la ligne de pantalon, en cm (variable \(X_6\)),
- de la hauteur de la tête aux pieds en cm (variable \(X_7\)),
- de la circonférence maximum du mollet en cm (variable \(X_8\)),
- de la circonférence de la cuisse, mesurée à mi-chemin entre le genou et le haut de la jambe, en cm (variable \(X_9\)),
- de la circonférence de la tête en cm (variable \(X_{10}\)).
On modélise le problème comme une rlm : \[ Y = \beta_0+\beta_1 X_1+\ldots+\beta_{10} X_{10} + \epsilon. \] où, \(\beta_0,\beta_1,\ldots,\beta_{10}\) sont 11 coefficients inconnus et \(\epsilon\sim\mathcal{N}(0,\sigma^2)\) avec \(\sigma\) inconnu.
Vous pouvez lire la section 6.5 Lab 1: Subset Selection Methods (6.5.1 and 6.5.2) page 244 du livre ISLR pour vous aider.
Modèle de régression linéaire multiple
1. Chargez les données PhyWeigth.
setwd("E:/CEPE DS/Session 2 regression")
PhyWeight <- read.table("PhyWeight.txt", header = T)2. Affichez les quatre premières lignes de la table.
head(PhyWeight,4)## Y X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
## 1 77.0 28.5 33.5 100 38.5 114 85.0 178 37.5 53 58
## 2 85.5 29.5 36.5 107 39.0 119 90.5 187 40.0 52 59
## 3 63.0 25.0 31.0 94 36.5 102 80.5 175 33.0 49 57
## 4 80.5 28.5 34.0 104 39.0 114 91.5 183 38.0 50 603. Qu’elle est la dimension de la table ?
dim(PhyWeight)## [1] 22 114. Affichez un résumé statistique des données.
str(PhyWeight)## 'data.frame': 22 obs. of 11 variables:
## $ Y : num 77 85.5 63 80.5 79.5 94 66 69 65 58 ...
## $ X1 : num 28.5 29.5 25 28.5 28.5 30.5 26.5 27 26.5 26.5 ...
## $ X2 : num 33.5 36.5 31 34 36.5 38 29 31 29 31 ...
## $ X3 : num 100 107 94 104 107 112 93 95 93 96 ...
## $ X4 : num 38.5 39 36.5 39 39 39 35 37 35 35 ...
## $ X5 : num 114 119 102 114 114 121 105 108 112 103 ...
## $ X6 : num 85 90.5 80.5 91.5 92 101 76 84 74 76 ...
## $ X7 : num 178 187 175 183 174 ...
## $ X8 : num 37.5 40 33 38 40 39.5 38.5 36 34 35 ...
## $ X9 : num 53 52 49 50 53 57.5 50 49 47 46 ...
## $ X10: num 58 59 57 60 59 59 58.5 60 55.5 58 ...summary(PhyWeight)## Y X1 X2 X3
## Min. :54.50 Min. :24.00 Min. :28.50 Min. : 87.50
## 1st Qu.:66.00 1st Qu.:26.50 1st Qu.:31.00 1st Qu.: 94.25
## Median :73.00 Median :28.25 Median :33.50 Median : 99.50
## Mean :73.93 Mean :27.77 Mean :33.20 Mean : 99.64
## 3rd Qu.:80.38 3rd Qu.:29.38 3rd Qu.:35.88 3rd Qu.:105.50
## Max. :94.00 Max. :31.00 Max. :38.00 Max. :112.00
## X4 X5 X6 X7
## Min. :35.00 Min. :102.0 Min. : 74.00 Min. :168.5
## 1st Qu.:35.75 1st Qu.:105.8 1st Qu.: 80.62 1st Qu.:173.2
## Median :38.50 Median :113.5 Median : 84.00 Median :178.8
## Mean :37.59 Mean :111.8 Mean : 85.57 Mean :178.5
## 3rd Qu.:39.00 3rd Qu.:117.8 3rd Qu.: 91.25 3rd Qu.:182.9
## Max. :40.50 Max. :121.0 Max. :101.00 Max. :188.0
## X8 X9 X10
## Min. :32.00 Min. :42.00 Min. :55.50
## 1st Qu.:36.00 1st Qu.:49.00 1st Qu.:57.00
## Median :38.00 Median :49.75 Median :58.00
## Mean :37.34 Mean :49.73 Mean :58.07
## 3rd Qu.:38.88 3rd Qu.:51.75 3rd Qu.:59.00
## Max. :42.00 Max. :57.50 Max. :60.005. La table contient-elle des données manquantes ?
sum(is.na(PhyWeight))## [1] 0is.na(PhyWeight)## Y X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
## [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [4,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [7,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [8,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [9,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [10,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [11,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [14,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [15,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [16,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [17,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [18,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [19,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [20,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [21,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [22,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEanyNA(PhyWeight)## [1] FALSE6. Sur la même fenêtre, représentez graphiquement chaque paire de variables.
plot(PhyWeight, cex = .5, col = "darkgrey")
7. Construisez un modèle linéaire à partir de l’ensemble des variables explicatives de PhyWeigth.
lm.mod <- lm(Y ~., data=PhyWeight)8. Produire un résumé statistique des résultats du modèle. Quelles sont les variables les plus significatives ? Peut-on envisager d’effectuer une sélection de variables ?
summary(lm.mod)##
## Call:
## lm(formula = Y ~ ., data = PhyWeight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5523 -0.9965 0.0461 1.0499 4.1719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -69.51714 29.03739 -2.394 0.035605 *
## X1 1.78182 0.85473 2.085 0.061204 .
## X2 0.15509 0.48530 0.320 0.755275
## X3 0.18914 0.22583 0.838 0.420132
## X4 -0.48184 0.72067 -0.669 0.517537
## X5 -0.02931 0.23943 -0.122 0.904769
## X6 0.66144 0.11648 5.679 0.000143 ***
## X7 0.31785 0.13037 2.438 0.032935 *
## X8 0.44589 0.41251 1.081 0.302865
## X9 0.29721 0.30510 0.974 0.350917
## X10 -0.91956 0.52009 -1.768 0.104735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.287 on 11 degrees of freedom
## Multiple R-squared: 0.9772, Adjusted R-squared: 0.9565
## F-statistic: 47.17 on 10 and 11 DF, p-value: 1.408e-079. Extraire du summary la proportion de variance expliquée.
attributes(summary(lm.mod))## $names
## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"
##
## $class
## [1] "summary.lm"summary(lm.mod)["r.squared"]## $r.squared
## [1] 0.9772107Comme la plus part des variables explicatives ne sont pas significatives, nous pouvons envisager d’effectuer une sélection parmi ces variables.
10. Donner des intervalles de confiance pour les estimations des coefficients.
confint(lm.mod)## 2.5 % 97.5 %
## (Intercept) -133.42800877 -5.6062615
## X1 -0.09943047 3.6630678
## X2 -0.91303791 1.2232187
## X3 -0.30791607 0.6861870
## X4 -2.06801804 1.1043439
## X5 -0.55628825 0.4976636
## X6 0.40506844 0.9178140
## X7 0.03090785 0.6047850
## X8 -0.46204518 1.3538255
## X9 -0.37430495 0.9687296
## X10 -2.06427506 0.2251497Best subset
11. Utilisez la fonction regsubsets() (du package leaps pour effectuer une sélection de variables via l’approche exhaustive Best Subset.
library(leaps)## Warning: package 'leaps' was built under R version 3.3.3regfit.full <- regsubsets(Y ~., PhyWeight, method="exhaustive", nvmax = 10)12. Utilisez la commande summary() pour afficher le meilleur ensemble de variables pour chaque taille de modèle (un astérisque indique qu’une variable donnée est incluse dans le modèle correspondant).
reg.summary <- summary(regfit.full)
reg.summary## Subset selection object
## Call: regsubsets.formula(Y ~ ., PhyWeight, method = "exhaustive", nvmax = 10)
## 10 Variables (and intercept)
## Forced in Forced out
## X1 FALSE FALSE
## X2 FALSE FALSE
## X3 FALSE FALSE
## X4 FALSE FALSE
## X5 FALSE FALSE
## X6 FALSE FALSE
## X7 FALSE FALSE
## X8 FALSE FALSE
## X9 FALSE FALSE
## X10 FALSE FALSE
## 1 subsets of each size up to 10
## Selection Algorithm: exhaustive
## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
## 1 ( 1 ) " " " " " " " " " " "*" " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " "*" " " " " " " " "
## 3 ( 1 ) "*" " " " " " " " " "*" "*" " " " " " "
## 4 ( 1 ) "*" " " " " " " " " "*" "*" " " "*" " "
## 5 ( 1 ) "*" " " " " " " " " "*" "*" " " "*" "*"
## 6 ( 1 ) "*" " " " " " " " " "*" "*" "*" "*" "*"
## 7 ( 1 ) "*" " " "*" " " " " "*" "*" "*" "*" "*"
## 8 ( 1 ) "*" " " "*" "*" " " "*" "*" "*" "*" "*"
## 9 ( 1 ) "*" "*" "*" "*" " " "*" "*" "*" "*" "*"
## 10 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"13. En considérant le \(R^2\), utilisez la fonction plot() pour tracer une table de modèles montrant quelles variables sont contenu dans chaque modèle (voir ?plot.regsubsets).
plot(regfit.full, scale="bic")
On constate que le plus grand \(R^2\) est obtenu avec la présence de toutes les variables car la valeur \(R^2\) d’un modèle augmente toujours lorsqu’un nouveau terme est ajouté. Un modèle contenant un plus grand nombre de termes peut, par ce simple fait, sembler présenter un meilleur ajustement. Le \(R^2\) ajusté permet de déterminer à quel point le modèle ajuste vos données lorsque vous souhaitez l’ajuster en fonction du nombre de prédicteurs inclus. La valeur du \(R^2\) ajusté intègre le nombre de prédicteurs dans le modèle elle donc plus adaptée pour nous aider à choisir le modèle.
14. Quelle est l’augmentation du \(R^2\) quand une seule variable est incluse dans le modèle ? Quand toutes les variables sont incluses ?
reg.summary$rsq## [1] 0.8373691 0.9364001 0.9555697 0.9659393 0.9707009 0.9743589 0.9762400
## [8] 0.9769473 0.9771796 0.9772107reg.summary$rsq[1]## [1] 0.8373691reg.summary$rsq[10]## [1] 0.977210715. Afin de nous aider à choisir le modèle à sélectionner, identifiez l’emplacement du point maximum / minimum pour chaque critère : \(RSS\), \(R^2\) ajusté, \(C_p\) et \(BIC\). Dans chaque cas, afficher les variables sélectionnées.
min.rss <- which.min(reg.summary$rss)
max.adjr2 <- which.max(reg.summary$adjr2)
min.cp <- which.min(reg.summary$cp)
min.bic <- which.min(reg.summary$bic)
min.rss## [1] 10max.adjr2## [1] 7min.cp## [1] 5min.bic## [1] 5names(which(reg.summary$which[min.rss,]==TRUE))## [1] "(Intercept)" "X1" "X2" "X3" "X4"
## [6] "X5" "X6" "X7" "X8" "X9"
## [11] "X10"names(which(reg.summary$which[max.adjr2,]==TRUE))## [1] "(Intercept)" "X1" "X3" "X6" "X7"
## [6] "X8" "X9" "X10"names(which(reg.summary$which[min.cp,]==TRUE))## [1] "(Intercept)" "X1" "X6" "X7" "X9"
## [6] "X10"names(which(reg.summary$which[min.bic,]==TRUE))## [1] "(Intercept)" "X1" "X6" "X7" "X9"
## [6] "X10"16. Sur une même fenêtre graphique représenter les courbes des différents critère. Ajouter sur chaque courbe, le maximum/minimum correspondant.
par(mfrow =c(2,2))
plot(reg.summary$rss,xlab="Number of Variables",ylab="RSS",type="l")
points(min.rss,reg.summary$rss[min.rss],col ="red",cex =2, pch =20)
plot(reg.summary$adjr2,xlab="Number of Variables ",ylab="Adjusted RSq",type="l")
points(max.adjr2,reg.summary$adjr2[max.adjr2],col ="red",cex =2, pch =20)
plot(reg.summary$cp,xlab="Number of Variables ",ylab="Cp",type="l")
points(min.cp,reg.summary$cp[min.cp],col ="red",cex =2, pch =20)
plot(reg.summary$bic,xlab="Number of Variables ",ylab="BIC",type="l")
points(min.bic,reg.summary$bic[min.bic],col ="red",cex =2, pch =20)
17. Affichez les variables sélectionnées pour le meilleur modèle avec un nombre donné de prédicteurs, classés selon chaque critère.
plot(regfit.full,scale ="r2")
plot(regfit.full,scale ="adjr2")
plot(regfit.full,scale ="Cp")
plot(regfit.full,scale ="bic")
18. Faire une régression avec le meilleur modèle selon la statistique BIC.
var.bic <- names(which(reg.summary$which[min.bic,]==TRUE))
var.bic.formula <- paste("Y", "~", paste(var.bic[-1], collapse=" + "))
var.bic.formula## [1] "Y ~ X1 + X6 + X7 + X9 + X10"best.model <- lm(var.bic.formula, data=PhyWeight)19. Qu’elle est la proportion de variance expliquée ? Comparer sa valeur à celle du modèle complet.
summary(best.model)##
## Call:
## lm(formula = var.bic.formula, data = PhyWeight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1865 -0.9735 -0.1115 0.9288 4.9834
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -80.45330 24.72024 -3.255 0.00497 **
## X1 2.12319 0.44541 4.767 0.00021 ***
## X6 0.66561 0.10040 6.630 5.79e-06 ***
## X7 0.27704 0.08179 3.387 0.00376 **
## X9 0.52317 0.22720 2.303 0.03506 *
## X10 -0.63714 0.39512 -1.613 0.12639
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.15 on 16 degrees of freedom
## Multiple R-squared: 0.9707, Adjusted R-squared: 0.9615
## F-statistic: 106 on 5 and 16 DF, p-value: 1.104e-11cat(paste0("La proportion de variance expliquée est de : ",
round(unlist(summary(best.model)["r.squared"]),2)))## La proportion de variance expliquée est de : 0.97cat(paste0("Celle du modèle complet est de : ",
round(unlist(summary(lm.mod)["r.squared"]),2)))## Celle du modèle complet est de : 0.9820. Afficher des informations sur le meilleur modèle pour vérifier les hypothèses standards du modèle linéaire (par exemple, l’hypothèse de normalité du bruit, l’indépendance du bruit, l’égalité de variance,…)
On trace les graphiques permettant de conclure sur la validité du modèle linéaire
par(mfrow = c(2, 2))
plot(best.model,1:4)
Les résultats obtenus sont consistents. En particulier :
le graphique en haut à gauche montre que le nuage de point semble répartit de manière aléatoire : le modèle linéaire semble bien adapté,
le graphique en haut à droite montre un alignement des points sur la diagonal ce qui nous permet d’admettre l’indépendance \(\epsilon_1,\ldots,\epsilon_n\),
celui en bas à gauche ne montre pas de structure particulière : on admet, l’égalité des variances de \(\epsilon_1,\ldots,\epsilon_n\),
celui en bas à droite montre les distances de Cook : aucune d’entres elles ne dépassent 1, il n’y a donc pas de valeurs anormales (cependant, la distance associée à l’individus 14 est élevée, on pourrait envisager de la supprimer de la table).
Forward et Backward Stepwise
21. Utilisez la fonction regsubsets() pour effectuer une sélection de variables en utilisant successivement les approches forward stepwise et backward stepwise.
regfit.fwd <- regsubsets(Y~.,data=PhyWeight, nvmax=10, method="forward")
reg.summary.fwd <- summary(regfit.fwd)
regfit.bwd <- regsubsets(Y~.,data=PhyWeight ,nvmax=10, method="backward")
reg.summary.bwd <- summary(regfit.bwd)22. Pour les deux méthodes, affichez les variables sélectionnées selon les critères \(C_p\) et \(BIC\). Les variables sélectionnées sont-elles identiques à celles sélectionnées par la recherche exhaustive ?
max.adjr2.fwd <- which.max(reg.summary.fwd$adjr2)
min.cp.fwd <- which.min(reg.summary.fwd$cp)
min.bic.fwd <- which.min(reg.summary.fwd$bic)
which(reg.summary.fwd$which[max.adjr2.fwd,]==T)## (Intercept) X1 X3 X6 X7 X8
## 1 2 4 7 8 9
## X9 X10
## 10 11which(reg.summary.fwd$which[min.cp.fwd,]==T)## (Intercept) X1 X6 X7 X9 X10
## 1 2 7 8 10 11which(reg.summary.fwd$which[min.bic.fwd,]==T)## (Intercept) X1 X6 X7 X9 X10
## 1 2 7 8 10 11max.adjr2.bwd <- which.max(reg.summary.bwd$adjr2)
min.cp.bwd <- which.min(reg.summary.bwd$cp)
min.bic.bwd <- which.min(reg.summary.bwd$bic)
which(reg.summary.bwd$which[max.adjr2.bwd,]==T)## (Intercept) X1 X3 X6 X7 X8
## 1 2 4 7 8 9
## X9 X10
## 10 11which(reg.summary.bwd$which[min.cp.bwd,]==T)## (Intercept) X1 X6 X7 X9 X10
## 1 2 7 8 10 11which(reg.summary.bwd$which[min.bic.bwd,]==T)## (Intercept) X1 X6 X7 X9 X10
## 1 2 7 8 10 11Approche apprentissage/test.
La performance du modèle issu d’une méthode d’apprentissage s’évalue par sa capacité de prévision. La mesure de cette performance est très importante puisque, d’une part, elle permet d’opérer une sélection de modèle dans une famille associée à la méthode d’apprentissage utilisée et, d’autre part, elle guide le choix de la méthode en comparant chacun des modèles optimisés à l’étape précédente. Enfin, elle fournit une mesure de la qualité ou encore de la confiance que l’on peut accorder à la prévision.
23. Partager l’échantillon en un ensemble d’apprentissage et un ensemble test (par exemple en prenant 2/3:1/3).
set.seed(10)
train <- sample(1:nrow(PhyWeight),2*nrow(PhyWeight)/3)
test <- (-train)
test## [1] -12 -7 -9 -14 -2 -4 -5 -16 -20 -6 -8 -21 -18 -1324. Utiliser regsubsets() sur l’ensemble d’apprentissage à l’aide de la méthode exhaustive.
regfit.best <- regsubsets(Y~.,data=PhyWeight[train,],nvmax=10)25. Calculez l’erreur de test pour le meilleur modèle de chaque taille.
test.mat <- model.matrix(Y~.,data=PhyWeight[test,])
# initialisation de l'erreur de prediction
val.errors <- rep(NA ,10)
for(i in 1:10){
# extraction des estimateurs des coefs
coefi <- coef(regfit.best,id=i)
# calcul de la prediction
pred <- test.mat[,names(coefi)]%*%coefi
# calcul de l'erreur de prediction
val.errors[i] <- mean((PhyWeight$Y[test]-pred)^2)
}
val.errors## [1] 107.22382 12.87755 23.56798 21.38280 19.25576 83.75554 114.39062
## [8] 140.27235 129.48728 145.68507which.min(val.errors)## [1] 2coef(regfit.best,which.min(val.errors))## (Intercept) X1 X6
## -64.278145 2.679237 0.75445126. Sur l’ensemble des données, effectuer une sélection de variables par la méthode exhaustive et selectionner le meilleur modèle.
regfit.best <- regsubsets(Y~.,data=PhyWeight,nvmax=10)
coef(regfit.best,which.min(val.errors))## (Intercept) X1 X6
## -68.7183553 2.7546175 0.773031927. Dans ce cas, est-ce que l’approche d’apprentissage/test est performante ? Expliquer brièvement pourquoi.
Non !
L’approche d’apprentissage/test possède quelques inconvénients :
l’estimation de validation du taux d’erreur de test peut être très variable (selon précisément quelles observations sont incluses dans la base d’apprentissage et quelles observations sont incluses dans le jeu de test),
seul un sous-ensemble des observations est utilisé pour ajuster le modèle.
un autre problème est causé par la petite taille de l’échantillon. Le nombre d’échantillons n’est généralement pas suffisant pour être représentatif de la population totale, ce qui est le cas ici.