TP Abidjan

SEANCE DE TP2

EXERCICE 1

Importons tout d’abord le jeu de données.

euc<-[Link]("~/Dropbox/ABIDJAN/DONNEES/[Link]",header=T)
head(euc)

## id ht circ
## 1 1 18.25 36
## 2 2 19.75 42
## 3 3 16.50 33
## 4 4 18.25 39
## 5 5 19.50 43
## 6 6 16.25 34
1. Nuage de points :
data<-euc[,-1] #On enlève la première colonne
plot(data$circ,data$ht)
25
data$ht

20
15

30 40 50 60 70

data$circ

Le nuage de points révèle un dépendance de type “racine carré”, on rajoute donc cette variable en tant que
variable explicative.

1
rac<-sqrt(data$circ)
donnees<-cbind(data,rac) #Concaténation des 2 tableaux

Construction des 3 modèles et summary. Par exemple, model1 correspond au modèle :

√
ht = β1 circ + β2 circ + ε,
avec ε ∼ N (0, σ 2 ), ce qui nous fait passer d’un modèle linéaire à un modèle non linéaire.

model1<-lm(ht~.,data=donnees)
model2<-lm(ht~circ,data=donnees)
model3<-lm(ht~rac,data=donnees)

summary(model1)

##
## Call:
## lm(formula = ht ~ ., data = donnees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1881 -0.6881 0.0427 0.7927 3.7481
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -24.35200 2.61444 -9.314 <2e-16 ***
## circ -0.48295 0.05793 -8.336 <2e-16 ***
## rac 9.98689 0.78033 12.798 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.136 on 1426 degrees of freedom
## Multiple R-squared: 0.7922, Adjusted R-squared: 0.7919
## F-statistic: 2718 on 2 and 1426 DF, p-value: < 2.2e-16
summary(model2)

##
## Call:
## lm(formula = ht ~ circ, data = donnees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7659 -0.7802 0.0557 0.8271 3.6913
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.037476 0.179802 50.26 <2e-16 ***
## circ 0.257138 0.003738 68.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.199 on 1427 degrees of freedom
## Multiple R-squared: 0.7683, Adjusted R-squared: 0.7682
## F-statistic: 4732 on 1 and 1427 DF, p-value: < 2.2e-16

2
summary(model3)

##
## Call:
## lm(formula = ht ~ rac, data = donnees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5360 -0.7249 0.0265 0.7813 3.6904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.73036 0.33600 -8.126 9.51e-16 ***
## rac 3.49424 0.04883 71.560 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.163 on 1427 degrees of freedom
## Multiple R-squared: 0.7821, Adjusted R-squared: 0.7819
## F-statistic: 5121 on 1 and 1427 DF, p-value: < 2.2e-16

Tests de nullité des coefficients supplémentaires dans les modèles emboîtés.

anova(model1,model2)

## Analysis of Variance Table

##
## Model 1: ht ~ circ + rac
## Model 2: ht ~ circ
## [Link] RSS Df Sum of Sq F Pr(>F)
## 1 1426 1840.7
## 2 1427 2052.1 -1 -211.43 163.8 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model1,model3)

## Analysis of Variance Table

##
## Model 1: ht ~ circ + rac
## Model 2: ht ~ rac
## [Link] RSS Df Sum of Sq F Pr(>F)
## 1 1426 1840.7
## 2 1427 1930.3 -1 -89.696 69.489 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Comparaison des 3 modèles par VC holdout (ne pas oublier de faire tout d’abord une permutation aléatoire
sur les lignes, au cas où les données auraient été rentrées dans un ordre) :

[Link](123)
n<-nrow(donnees)
perm<-sample(n)

3
donnees<-donnees[perm[1:n],] #Permutation aléatoire sur les lignes

dapp<-donnees[1:1200,] #Construction des échantillons d'apprentissage et test

dtest<-donnees[-(1:1200),]

model1_app<-lm(ht~.,data=dapp) #Construction des 3 modèles sur l'échantillon d'apprentissage

model2_app<-lm(ht~circ,data=dapp)
model3_app<-lm(ht~rac,data=dapp)

pred1<-predict(model1_app,newdata=dtest,type="response") #Calculs des prédictions sur l'échantillon t

pred2<-predict(model2_app,newdata=dtest,type="response")
pred3<-predict(model3_app,newdata=dtest,type="response")

mean((dtest$ht-pred1)ˆ2) #Calculs des erreurs de prédictions sur l'échantillon test

## [1] 1.283919
mean((dtest$ht-pred2)ˆ2)

## [1] 1.434178
mean((dtest$ht-pred3)ˆ2)

## [1] 1.349997

3. Mais pour plus de fiabilité du résultat, il faut refaire cette procédure plusieurs fois. On peut aussi envisager
une validation croisée K-fold pour les 3 types de modèles :

library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 4.1-3
library(boot)

#Construction des modèles : pour le VC K-fold, il faut utiliser la fonction glm (Generalized Linear Mode
model1_bis<-glm(ht~.,donnees,family="gaussian")
model2_bis<-glm(ht~circ,donnees,family="gaussian")
model3_bis<-glm(ht~rac,donnees,family="gaussian")

[Link](donnees,model1_bis,K=10)$delta[1] #Calculs des erreurs quadratiques par VC 10-fold

## [1] 1.294896
[Link](donnees,model2_bis,K=10)$delta[1]

## [1] 1.438287
[Link](donnees,model3_bis,K=10)$delta[1]

## [1] 1.353946

On retient donc model1.

4. Il reste à vérifier que ce modèle est valide en examinant ses résidus.

4
plot(model1_bis)

4 Residuals vs Fitted

15 1129
2
Residuals

0
−2
−4

49

14 16 18 20 22 24 26

Predicted values
glm(ht ~ .)
Normal Q−Q
4

1129 15
Std. Pearson resid.

2
0
−2

49

−3 −2 −1 0 1 2 3

Theoretical Quantiles
glm(ht ~ .)

5
 Scale−Location
49
15 1129
1.5
Std. Pearson resid.

1.0
0.5
0.0

14 16 18 20 22 24 26

Predicted values
glm(ht ~ .)
Residuals vs Leverage
4

15
2
Std. Pearson resid.

1233
0
−2

318
−4

Cook's distance

0.00 0.01 0.02 0.03 0.04

Leverage
glm(ht ~ .)

Eventuellement, pour être plus précis, on peut lancer un test de normalité, un test de décorrélation des
résidus, un test de linéarité et enfin un test d’homogénéité de la variance des résidus.

library(zoo)

##

6
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## [Link], [Link]
library(lmtest)
raintest(model1) #Test de linéarité du modèle

##
## Rainbow test
##
## data: model1
## Rain = 1.1144, df1 = 715, df2 = 711, p-value = 0.07423
[Link](resid(model1)) #Test de normalité des résidus

##
## Shapiro-Wilk normality test
##
## data: resid(model1)
## W = 0.99802, p-value = 0.08508
gqtest(model1) #Test d'homogénéité de la variance des résidus

##
## Goldfeld-Quandt test
##
## data: model1
## GQ = 1.1471, df1 = 712, df2 = 711, p-value = 0.03372
## alternative hypothesis: variance increases from segment 1 to 2
dwtest(model1) #Test de décorrélation des résidus

##
## Durbin-Watson test
##
## data: model1
## DW = 0.94725, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

5. On sait que la meilleure prévision possible, au sens du critère quadratique, est en notant β̂1 et β̂2 les
estimations des coefficients :

√
ht = β̂0 + β̂1 circ + β̂2 circ,
et la variable aléatoire ε sert alors à ajuster les bornes de l’intervalle de confiance.

On prédit la hauteur d’un nouvel arbre avec circ=50. On prend model1 issu de la fonction lm :

predict(model1,[Link](rac=sqrt(50),circ=50), interval="prediction",level=0.95)

## fit lwr upr

## 1 22.11869 19.88891 24.34847

Mais avec model1bis qui est de la classe “glm”, on n’a que des valeurs ajustées :

7
predict(model1_bis,newdata=[Link](rac=sqrt(50),circ=50))

## 1
## 22.11869

8
SEANCE DE TP3
EXERCICE 1

La base de données comprend les mesures sur 50 états des USA de la population en 1975, du revenu par
individu (income) en 1974, du pourcentage d’illetrisme (illiteracy) en 1970, de l’espérance de vie moyennes
([Link]) sur la période 1969-1971,. . .

state<-[Link](state.x77,[Link]=[Link])
head(state)

## Population Income Illiteracy [Link] Murder [Link] Frost Area

## AL 3615 3624 2.1 69.05 15.1 41.3 20 50708
## AK 365 6315 1.5 69.31 11.3 66.7 152 566432
## AZ 2212 4530 1.8 70.55 7.8 58.1 15 113417
## AR 2110 3378 1.9 70.66 10.1 39.9 65 51945
## CA 21198 5114 1.1 71.71 10.3 62.6 20 156361
## CO 2541 4884 0.7 72.06 6.8 63.9 166 103766

Avant toute chose, puisque nous allons effectuer dans la question 2 une recherche exhaustive du meilleur
modèle, il faut placer [Link] en dernière colonne :

state<-state[,c(1:3,5:8,4)]
head(state)

## Population Income Illiteracy Murder [Link] Frost Area [Link]

## AL 3615 3624 2.1 15.1 41.3 20 50708 69.05
## AK 365 6315 1.5 11.3 66.7 152 566432 69.31
## AZ 2212 4530 1.8 7.8 58.1 15 113417 70.55
## AR 2110 3378 1.9 10.1 39.9 65 51945 70.66
## CA 21198 5114 1.1 10.3 62.6 20 156361 71.71
## CO 2541 4884 0.7 6.8 63.9 166 103766 72.06

1. Calculons tout d’abord le modèle complet.

[Link]<-lm([Link]~.,data=state)
summary([Link])

##
## Call:
## lm(formula = [Link] ~ ., data = state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.48895 -0.51232 -0.02747 0.57002 1.49447
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.094e+01 1.748e+00 40.586 < 2e-16 ***
## Population 5.180e-05 2.919e-05 1.775 0.0832 .
## Income -2.180e-05 2.444e-04 -0.089 0.9293
## Illiteracy 3.382e-02 3.663e-01 0.092 0.9269
## Murder -3.011e-01 4.662e-02 -6.459 8.68e-08 ***

9
## [Link] 4.893e-02 2.332e-02 2.098 0.0420 *
## Frost -5.735e-03 3.143e-03 -1.825 0.0752 .
## Area -7.383e-08 1.668e-06 -0.044 0.9649
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7448 on 42 degrees of freedom
## Multiple R-squared: 0.7362, Adjusted R-squared: 0.6922
## F-statistic: 16.74 on 7 and 42 DF, p-value: 2.534e-10

2. Appliquons une méthode de sélection exhaustive, selon AIC et BIC :

library(leaps)
library(bestglm)
[Link]<-bestglm(state,family=gaussian,IC="AIC")
[Link]<-bestglm(state,family=gaussian,IC="BIC")
summary([Link]$BestModel)

##
## Call:
## lm(formula = y ~ ., data = [Link](Xy[, c(bestset[-1], FALSE),
## drop = FALSE], y = y))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## Population 5.014e-05 2.512e-05 1.996 0.05201 .
## Murder -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## [Link] 4.658e-02 1.483e-02 3.142 0.00297 **
## Frost -5.943e-03 2.421e-03 -2.455 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12
summary([Link]$BestModel)

##
## Call:
## lm(formula = y ~ ., data = [Link](Xy[, c(bestset[-1], FALSE),
## drop = FALSE], y = y))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)

10
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## Population 5.014e-05 2.512e-05 1.996 0.05201 .
## Murder -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## [Link] 4.658e-02 1.483e-02 3.142 0.00297 **
## Frost -5.943e-03 2.421e-03 -2.455 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12

3. Calculons maintenant le meilleur modèle obtenu par les méthodes pas à pas descendantes. Par défaut, le
critère AIC :

[Link]<-step([Link],data=state,direction="backward")

## Start: AIC=-22.18
## [Link] ~ Population + Income + Illiteracy + Murder + [Link] +
## Frost + Area
##
## Df Sum of Sq RSS AIC
## - Area 1 0.0011 23.298 -24.182
## - Income 1 0.0044 23.302 -24.175
## - Illiteracy 1 0.0047 23.302 -24.174
## <none> 23.297 -22.185
## - Population 1 1.7472 25.044 -20.569
## - Frost 1 1.8466 25.144 -20.371
## - [Link] 1 2.4413 25.738 -19.202
## - Murder 1 23.1411 46.438 10.305
##
## Step: AIC=-24.18
## [Link] ~ Population + Income + Illiteracy + Murder + [Link] +
## Frost
##
## Df Sum of Sq RSS AIC
## - Illiteracy 1 0.0038 23.302 -26.174
## - Income 1 0.0059 23.304 -26.170
## <none> 23.298 -24.182
## - Population 1 1.7599 25.058 -22.541
## - Frost 1 2.0488 25.347 -21.968
## - [Link] 1 2.9804 26.279 -20.163
## - Murder 1 26.2721 49.570 11.569
##
## Step: AIC=-26.17
## [Link] ~ Population + Income + Murder + [Link] + Frost
##
## Df Sum of Sq RSS AIC
## - Income 1 0.006 23.308 -28.161
## <none> 23.302 -26.174
## - Population 1 1.887 25.189 -24.280
## - Frost 1 3.037 26.339 -22.048
## - [Link] 1 3.495 26.797 -21.187
## - Murder 1 34.739 58.041 17.456

11
##
## Step: AIC=-28.16
## [Link] ~ Population + Murder + [Link] + Frost
##
## Df Sum of Sq RSS AIC
## <none> 23.308 -28.161
## - Population 1 2.064 25.372 -25.920
## - Frost 1 3.122 26.430 -23.877
## - [Link] 1 5.112 28.420 -20.246
## - Murder 1 34.816 58.124 15.528
summary([Link])

##
## Call:
## lm(formula = [Link] ~ Population + Murder + [Link] + Frost,
## data = state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## Population 5.014e-05 2.512e-05 1.996 0.05201 .
## Murder -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## [Link] 4.658e-02 1.483e-02 3.142 0.00297 **
## Frost -5.943e-03 2.421e-03 -2.455 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12

Pour le critère BIC, il faut préciser la valeur de la pénalité “k=log(nrow(state))” (noter que les sorties notent
encore “AIC”, alors que c’est bien le BIC qui est utilisé).

[Link]<-step([Link],direction="backward",k=log(nrow(state)))

## Start: AIC=-6.89
## [Link] ~ Population + Income + Illiteracy + Murder + [Link] +
## Frost + Area
##
## Df Sum of Sq RSS AIC
## - Area 1 0.0011 23.298 -10.7981
## - Income 1 0.0044 23.302 -10.7910
## - Illiteracy 1 0.0047 23.302 -10.7903
## - Population 1 1.7472 25.044 -7.1846
## - Frost 1 1.8466 25.144 -6.9866
## <none> 23.297 -6.8884
## - [Link] 1 2.4413 25.738 -5.8178
## - Murder 1 23.1411 46.438 23.6891

12
##
## Step: AIC=-10.8
## [Link] ~ Population + Income + Illiteracy + Murder + [Link] +
## Frost
##
## Df Sum of Sq RSS AIC
## - Illiteracy 1 0.0038 23.302 -14.7021
## - Income 1 0.0059 23.304 -14.6975
## - Population 1 1.7599 25.058 -11.0691
## <none> 23.298 -10.7981
## - Frost 1 2.0488 25.347 -10.4960
## - [Link] 1 2.9804 26.279 -8.6912
## - Murder 1 26.2721 49.570 23.0406
##
## Step: AIC=-14.7
## [Link] ~ Population + Income + Murder + [Link] + Frost
##
## Df Sum of Sq RSS AIC
## - Income 1 0.006 23.308 -18.601
## - Population 1 1.887 25.189 -14.720
## <none> 23.302 -14.702
## - Frost 1 3.037 26.339 -12.488
## - [Link] 1 3.495 26.797 -11.627
## - Murder 1 34.739 58.041 27.017
##
## Step: AIC=-18.6
## [Link] ~ Population + Murder + [Link] + Frost
##
## Df Sum of Sq RSS AIC
## <none> 23.308 -18.601
## - Population 1 2.064 25.372 -18.271
## - Frost 1 3.122 26.430 -16.228
## - [Link] 1 5.112 28.420 -12.598
## - Murder 1 34.816 58.124 23.176
summary([Link])

##
## Call:
## lm(formula = [Link] ~ Population + Murder + [Link] + Frost,
## data = state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## Population 5.014e-05 2.512e-05 1.996 0.05201 .
## Murder -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## [Link] 4.658e-02 1.483e-02 3.142 0.00297 **
## Frost -5.943e-03 2.421e-03 -2.455 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12

Pour des procédures ascendantes (avec l’AIC), il faut spécifier le scope, depuis le modèle initial (e.g. reg0)
vers le modèle final (e.g. [Link]):

reg0<-lm([Link]~1.,data=state)
[Link]<-step(reg0,scope=formula([Link]),direction="forward")

## Start: AIC=30.44
## [Link] ~ 1
##
## Df Sum of Sq RSS AIC
## + Murder 1 53.838 34.461 -14.609
## + Illiteracy 1 30.578 57.721 11.179
## + [Link] 1 29.931 58.368 11.737
## + Income 1 10.223 78.076 26.283
## + Frost 1 6.064 82.235 28.878
## <none> 88.299 30.435
## + Area 1 1.017 87.282 31.856
## + Population 1 0.409 87.890 32.203
##
## Step: AIC=-14.61
## [Link] ~ Murder
##
## Df Sum of Sq RSS AIC
## + [Link] 1 4.6910 29.770 -19.925
## + Population 1 4.0161 30.445 -18.805
## + Frost 1 3.1346 31.327 -17.378
## + Income 1 2.4047 32.057 -16.226
## <none> 34.461 -14.609
## + Area 1 0.4697 33.992 -13.295
## + Illiteracy 1 0.2732 34.188 -13.007
##
## Step: AIC=-19.93
## [Link] ~ Murder + [Link]
##
## Df Sum of Sq RSS AIC
## + Frost 1 4.3987 25.372 -25.920
## + Population 1 3.3405 26.430 -23.877
## <none> 29.770 -19.925
## + Illiteracy 1 0.4419 29.328 -18.673
## + Area 1 0.2775 29.493 -18.394
## + Income 1 0.1022 29.668 -18.097
##
## Step: AIC=-25.92
## [Link] ~ Murder + [Link] + Frost
##
## Df Sum of Sq RSS AIC
## + Population 1 2.06358 23.308 -28.161
## <none> 25.372 -25.920
## + Income 1 0.18232 25.189 -24.280

14
## + Illiteracy 1 0.17184 25.200 -24.259
## + Area 1 0.02573 25.346 -23.970
##
## Step: AIC=-28.16
## [Link] ~ Murder + [Link] + Frost + Population
##
## Df Sum of Sq RSS AIC
## <none> 23.308 -28.161
## + Income 1 0.0060582 23.302 -26.174
## + Illiteracy 1 0.0039221 23.304 -26.170
## + Area 1 0.0007900 23.307 -26.163
summary([Link])

##
## Call:
## lm(formula = [Link] ~ Murder + [Link] + Frost + Population,
## data = state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## Murder -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## [Link] 4.658e-02 1.483e-02 3.142 0.00297 **
## Frost -5.943e-03 2.421e-03 -2.455 0.01802 *
## Population 5.014e-05 2.512e-05 1.996 0.05201 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12

Idem avec le BIC :

[Link]<-step(reg0,scope=formula([Link]),direction="forward",k=log(nrow(state)))

## Start: AIC=32.35
## [Link] ~ 1
##
## Df Sum of Sq RSS AIC
## + Murder 1 53.838 34.461 -10.785
## + Illiteracy 1 30.578 57.721 15.004
## + [Link] 1 29.931 58.368 15.561
## + Income 1 10.223 78.076 30.107
## <none> 88.299 32.347
## + Frost 1 6.064 82.235 32.702
## + Area 1 1.017 87.282 35.680
## + Population 1 0.409 87.890 36.027
##

15
## Step: AIC=-10.79
## [Link] ~ Murder
##
## Df Sum of Sq RSS AIC
## + [Link] 1 4.6910 29.770 -14.1894
## + Population 1 4.0161 30.445 -13.0687
## + Frost 1 3.1346 31.327 -11.6415
## <none> 34.461 -10.7852
## + Income 1 2.4047 32.057 -10.4900
## + Area 1 0.4697 33.992 -7.5593
## + Illiteracy 1 0.2732 34.188 -7.2712
##
## Step: AIC=-14.19
## [Link] ~ Murder + [Link]
##
## Df Sum of Sq RSS AIC
## + Frost 1 4.3987 25.372 -18.271
## + Population 1 3.3405 26.430 -16.228
## <none> 29.770 -14.189
## + Illiteracy 1 0.4419 29.328 -11.025
## + Area 1 0.2775 29.493 -10.746
## + Income 1 0.1022 29.668 -10.449
##
## Step: AIC=-18.27
## [Link] ~ Murder + [Link] + Frost
##
## Df Sum of Sq RSS AIC
## + Population 1 2.06358 23.308 -18.601
## <none> 25.372 -18.271
## + Income 1 0.18232 25.189 -14.720
## + Illiteracy 1 0.17184 25.200 -14.699
## + Area 1 0.02573 25.346 -14.410
##
## Step: AIC=-18.6
## [Link] ~ Murder + [Link] + Frost + Population
##
## Df Sum of Sq RSS AIC
## <none> 23.308 -18.601
## + Income 1 0.0060582 23.302 -14.702
## + Illiteracy 1 0.0039221 23.304 -14.697
## + Area 1 0.0007900 23.307 -14.691
summary([Link])

##
## Call:
## lm(formula = [Link] ~ Murder + [Link] + Frost + Population,
## data = state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)

16
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## Murder -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## [Link] 4.658e-02 1.483e-02 3.142 0.00297 **
## Frost -5.943e-03 2.421e-03 -2.455 0.01802 *
## Population 5.014e-05 2.512e-05 1.996 0.05201 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12

4. Les variables qui expliquent [Link] sont, pour tous les critères et les méthodes : Murder, [Link], Frost
et Population. Le modèle retenu est donc

[Link]<-lm([Link]~Murder+[Link]+Frost+Population,data=state)

5. Vérifions maintenant que les résidus du modèle [Link] sont conformes aux hypothèses du modèle linéaire
gaussien.

plot([Link])

Residuals vs Fitted
1.5

HI
0.5
Residuals

−0.5

SC
−1.5

ME

69 70 71 72

Fitted values
lm([Link] ~ Murder + [Link] + Frost + Population)

17
 Normal Q−Q

2 HI
Standardized residuals

1
0
−1

SC
−2

ME

−2 −1 0 1 2

Theoretical Quantiles
lm([Link] ~ Murder + [Link] + Frost + Population)
Scale−Location
HI
1.5

ME
Standardized residuals

SC
1.0
0.5
0.0

69 70 71 72

Fitted values
lm([Link] ~ Murder + [Link] + Frost + Population)

18
 Residuals vs Leverage
1
2 HI
0.5
Standardized residuals

1
0
−1

WA
−2

0.5
Cook'sME
distance

0.0 0.1 0.2 0.3 0.4

Leverage
lm([Link] ~ Murder + [Link] + Frost + Population)
raintest([Link]) #test de linéarité

##
## Rainbow test
##
## data: [Link]
## Rain = 1.1324, df1 = 25, df2 = 20, p-value = 0.3924
[Link](resid([Link])) #test de normalité

##
## Shapiro-Wilk normality test
##
## data: resid([Link])
## W = 0.97935, p-value = 0.525
gqtest([Link]) #test d'homogénéité de la variance des résidus

##
## Goldfeld-Quandt test
##
## data: [Link]
## GQ = 0.71805, df1 = 20, df2 = 20, p-value = 0.7672
## alternative hypothesis: variance increases from segment 1 to 2
dwtest([Link]) #test d'indépendance des résidus

##
## Durbin-Watson test
##
## data: [Link]
## DW = 1.914, p-value = 0.3854

19
## alternative hypothesis: true autocorrelation is greater than 0

Cette étude des résidus valide l’utilisation de [Link].

Remarque : on pourrait être tenté, dans un premier temps, de comparer les performances respectives de
[Link] et [Link] (par exemple) en calculant les erreurs de prévision. Il faudrait pour cela découper
l’échantillon en un échantillon d’apprentissage et un échantillon test. Or, chaque individu de l’échantillon est
un état des USA : le modèle construit avec l’échantillon d’apprentissage ne contiendrait que certains états, et
on évaluerait ses performances avec l’échantillon test sur d’autres états ! Cette méthode est donc à proscrire.

6. Pour finir, calculons l’espérance de vie pour un nouveau pays tel que population=4250, frost=80,
[Link]=75 et Murder=8 :

predict([Link],[Link](Murder=8,[Link]=75,Frost=80,Population=4250),
interval="prediction",level=0.95)

## fit lwr upr

## 1 71.85724 70.23526 73.47921

20
SEANCE DE TP4

EXERCICE 1

bodyfat <- [Link]("~/Dropbox/ABIDJAN/DONNEES/[Link]", header=T)

head(bodyfat)

## Density [Link] Age Weight Height Neck Chest Abdomen Waist Hip Thigh
## 237 1.0705 12.4 64 155.25 69.5 37.9 95.8 82.8 32.59842 94.5 61.2
## 175 1.0689 13.1 37 151.00 67.0 35.3 92.6 83.2 32.75591 96.4 60.0
## 35 1.0263 32.3 41 247.25 73.5 42.1 117.0 115.6 45.51181 116.1 71.2
## 45 1.0742 10.8 40 133.50 67.5 33.6 88.2 73.7 29.01575 88.5 53.3
## 39 1.0258 32.6 50 203.00 67.0 40.2 114.8 108.1 42.55906 102.5 61.3
## 170 1.0983 0.7 35 125.75 65.5 34.0 90.8 75.0 29.52756 89.2 50.0
## Knee Ankle Bicep Forearm Wrist
## 237 39.1 22.3 29.8 28.9 18.3
## 175 38.1 22.0 31.5 26.6 16.7
## 35 43.3 26.3 37.3 31.7 19.7
## 45 34.5 22.5 27.9 26.2 17.3
## 39 41.1 24.7 34.1 31.0 18.3
## 170 34.8 22.0 24.8 25.9 16.9

1. Ecriture du modèle complet :

[Link]<-glm(Weight~.,data=bodyfat,family="gaussian")
summary([Link])

##
## Call:
## glm(formula = Weight ~ ., family = "gaussian", data = bodyfat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.0780 -2.2159 -0.4946 1.7531 9.5216
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.394e+02 1.268e+02 -2.677 0.0113 *
## Density 1.042e+01 1.242e+02 0.084 0.9336
## [Link] 2.707e-02 3.057e-01 0.089 0.9300
## Age -1.476e-01 6.238e-02 -2.367 0.0238 *
## Height 2.316e+00 3.409e-01 6.793 8.20e-08 ***
## Neck 9.887e-02 5.636e-01 0.175 0.8618
## Chest 1.190e+00 2.009e-01 5.922 1.09e-06 ***
## Abdomen 6.197e+05 8.323e+05 0.745 0.4616
## Waist -1.574e+06 2.114e+06 -0.745 0.4616
## Hip 7.651e-01 2.858e-01 2.677 0.0114 *
## Thigh 4.191e-01 2.965e-01 1.413 0.1666
## Knee -7.157e-01 6.150e-01 -1.164 0.2527
## Ankle 1.076e+00 7.092e-01 1.517 0.1384
## Bicep 6.836e-01 4.121e-01 1.659 0.1063
## Forearm 9.051e-01 8.604e-01 1.052 0.3002
## Wrist 1.690e+00 1.166e+00 1.449 0.1564

21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.35473)
##
## Null deviance: 38290.40 on 49 degrees of freedom
## Residual deviance: 420.06 on 34 degrees of freedom
## AIC: 282.31
##
## Number of Fisher Scoring iterations: 2

Puis un modèle sélectionné avec la méthode pas à pas descendante et le critère AIC :

[Link]<-step([Link],direction="backward")

## Start: AIC=282.31
## Weight ~ Density + [Link] + Age + Height + Neck + Chest + Abdomen +
## Waist + Hip + Thigh + Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Density 1 420.15 280.32
## - [Link] 1 420.16 280.32
## - Neck 1 420.44 280.36
## - Waist 1 426.91 281.12
## - Abdomen 1 426.91 281.12
## - Forearm 1 433.73 281.91
## - Knee 1 436.79 282.27
## <none> 420.06 282.31
## - Thigh 1 444.74 283.17
## - Wrist 1 446.01 283.31
## - Ankle 1 448.50 283.59
## - Bicep 1 454.06 284.20
## - Age 1 489.26 287.94
## - Hip 1 508.59 289.88
## - Chest 1 853.28 315.75
## - Height 1 990.12 323.18
##
## Step: AIC=280.32
## Weight ~ [Link] + Age + Height + Neck + Chest + Abdomen + Waist +
## Hip + Thigh + Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - [Link] 1 420.16 278.32
## - Neck 1 420.48 278.36
## - Waist 1 426.92 279.12
## - Abdomen 1 426.92 279.12
## - Forearm 1 433.95 279.94
## - Knee 1 437.06 280.30
## <none> 420.15 280.32
## - Thigh 1 444.85 281.18
## - Wrist 1 449.93 281.75
## - Ankle 1 453.85 282.18
## - Bicep 1 454.93 282.30
## - Age 1 489.28 285.94

22
## - Hip 1 510.76 288.09
## - Chest 1 926.40 317.86
## - Height 1 1064.28 324.80
##
## Step: AIC=278.32
## Weight ~ Age + Height + Neck + Chest + Abdomen + Waist + Hip +
## Thigh + Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Neck 1 420.54 276.37
## - Waist 1 426.93 277.12
## - Abdomen 1 426.93 277.12
## - Forearm 1 436.12 278.19
## - Knee 1 437.08 278.30
## <none> 420.16 278.32
## - Thigh 1 448.40 279.58
## - Wrist 1 450.08 279.76
## - Ankle 1 455.15 280.32
## - Bicep 1 456.07 280.43
## - Age 1 489.31 283.94
## - Hip 1 518.37 286.83
## - Chest 1 926.75 315.88
## - Height 1 1113.31 325.05
##
## Step: AIC=276.37
## Weight ~ Age + Height + Chest + Abdomen + Waist + Hip + Thigh +
## Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Waist 1 427.19 275.15
## - Abdomen 1 427.19 275.15
## - Knee 1 437.58 276.36
## <none> 420.54 276.37
## - Forearm 1 443.74 277.05
## - Thigh 1 448.59 277.60
## - Wrist 1 451.50 277.92
## - Ankle 1 455.16 278.32
## - Bicep 1 456.75 278.50
## - Age 1 492.12 282.23
## - Hip 1 518.70 284.86
## - Chest 1 938.47 314.51
## - Height 1 1178.30 325.88
##
## Step: AIC=275.15
## Weight ~ Age + Height + Chest + Abdomen + Hip + Thigh + Knee +
## Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Knee 1 444.52 275.14
## <none> 427.19 275.15
## - Thigh 1 454.94 276.30
## - Bicep 1 457.79 276.61
## - Ankle 1 458.79 276.72
## - Forearm 1 461.39 277.01

23
## - Wrist 1 462.71 277.15
## - Age 1 504.45 281.47
## - Hip 1 534.06 284.32
## - Abdomen 1 553.04 286.06
## - Chest 1 938.63 312.51
## - Height 1 1181.75 324.03
##
## Step: AIC=275.14
## Weight ~ Age + Height + Chest + Abdomen + Hip + Thigh + Ankle +
## Bicep + Forearm + Wrist
##
## Df Deviance AIC
## <none> 444.52 275.14
## - Thigh 1 464.76 275.37
## - Ankle 1 465.05 275.40
## - Bicep 1 470.46 275.98
## - Wrist 1 478.54 276.83
## - Forearm 1 479.43 276.92
## - Hip 1 539.66 282.84
## - Age 1 545.31 283.36
## - Abdomen 1 578.32 286.30
## - Chest 1 940.40 310.61
## - Height 1 1189.98 322.38

Les variables pertinentes sont donc, de ce pont de vue : Age + Height + Chest + Abdomen + Hip + Thigh
+ Ankle + Bicep + Forearm + Wrist

Par une méthode pas à pas descendante avec le BIC :

[Link]<-step([Link],direction="backward",k=log(nrow(bodyfat)))

## Start: AIC=312.91
## Weight ~ Density + [Link] + Age + Height + Neck + Chest + Abdomen +
## Waist + Hip + Thigh + Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Density 1 420.15 309.00
## - [Link] 1 420.16 309.00
## - Neck 1 420.44 309.04
## - Waist 1 426.91 309.80
## - Abdomen 1 426.91 309.80
## - Forearm 1 433.73 310.59
## - Knee 1 436.79 310.95
## - Thigh 1 444.74 311.85
## - Wrist 1 446.01 311.99
## - Ankle 1 448.50 312.27
## - Bicep 1 454.06 312.88
## <none> 420.06 312.91
## - Age 1 489.26 316.62
## - Hip 1 508.59 318.56
## - Chest 1 853.28 344.43
## - Height 1 990.12 351.86
##

24
## Step: AIC=309
## Weight ~ [Link] + Age + Height + Neck + Chest + Abdomen + Waist +
## Hip + Thigh + Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - [Link] 1 420.16 305.09
## - Neck 1 420.48 305.13
## - Waist 1 426.92 305.89
## - Abdomen 1 426.92 305.89
## - Forearm 1 433.95 306.71
## - Knee 1 437.06 307.06
## - Thigh 1 444.85 307.95
## - Wrist 1 449.93 308.52
## - Ankle 1 453.85 308.95
## <none> 420.15 309.00
## - Bicep 1 454.93 309.07
## - Age 1 489.28 312.71
## - Hip 1 510.76 314.86
## - Chest 1 926.40 344.63
## - Height 1 1064.28 351.56
##
## Step: AIC=305.09
## Weight ~ Age + Height + Neck + Chest + Abdomen + Waist + Hip +
## Thigh + Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Neck 1 420.54 301.23
## - Waist 1 426.93 301.98
## - Abdomen 1 426.93 301.98
## - Forearm 1 436.12 303.05
## - Knee 1 437.08 303.15
## - Thigh 1 448.40 304.43
## - Wrist 1 450.08 304.62
## <none> 420.16 305.09
## - Ankle 1 455.15 305.18
## - Bicep 1 456.07 305.28
## - Age 1 489.31 308.80
## - Hip 1 518.37 311.68
## - Chest 1 926.75 340.73
## - Height 1 1113.31 349.90
##
## Step: AIC=301.23
## Weight ~ Age + Height + Chest + Abdomen + Waist + Hip + Thigh +
## Knee + Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Waist 1 427.19 298.10
## - Abdomen 1 427.19 298.10
## - Knee 1 437.58 299.30
## - Forearm 1 443.74 300.00
## - Thigh 1 448.59 300.54
## - Wrist 1 451.50 300.87
## <none> 420.54 301.23
## - Ankle 1 455.16 301.27

25
## - Bicep 1 456.75 301.44
## - Age 1 492.12 305.17
## - Hip 1 518.70 307.80
## - Chest 1 938.47 337.45
## - Height 1 1178.30 348.83
##
## Step: AIC=298.1
## Weight ~ Age + Height + Chest + Abdomen + Hip + Thigh + Knee +
## Ankle + Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Knee 1 444.52 296.17
## - Thigh 1 454.94 297.33
## - Bicep 1 457.79 297.65
## - Ankle 1 458.79 297.75
## - Forearm 1 461.39 298.04
## <none> 427.19 298.10
## - Wrist 1 462.71 298.18
## - Age 1 504.45 302.50
## - Hip 1 534.06 305.35
## - Abdomen 1 553.04 307.10
## - Chest 1 938.63 333.55
## - Height 1 1181.75 345.06
##
## Step: AIC=296.17
## Weight ~ Age + Height + Chest + Abdomen + Hip + Thigh + Ankle +
## Bicep + Forearm + Wrist
##
## Df Deviance AIC
## - Thigh 1 464.76 294.49
## - Ankle 1 465.05 294.52
## - Bicep 1 470.46 295.10
## - Wrist 1 478.54 295.95
## - Forearm 1 479.43 296.04
## <none> 444.52 296.17
## - Hip 1 539.66 301.96
## - Age 1 545.31 302.48
## - Abdomen 1 578.32 305.42
## - Chest 1 940.40 329.73
## - Height 1 1189.98 341.50
##
## Step: AIC=294.49
## Weight ~ Age + Height + Chest + Abdomen + Hip + Ankle + Bicep +
## Forearm + Wrist
##
## Df Deviance AIC
## - Wrist 1 487.36 292.95
## - Ankle 1 495.27 293.76
## - Bicep 1 499.94 294.23
## <none> 464.76 294.49
## - Forearm 1 524.89 296.66
## - Age 1 582.79 301.89
## - Abdomen 1 604.94 303.76
## - Hip 1 661.53 308.23

26
## - Chest 1 945.78 326.10
## - Height 1 1190.69 337.62
##
## Step: AIC=292.95
## Weight ~ Age + Height + Chest + Abdomen + Hip + Ankle + Bicep +
## Forearm
##
## Df Deviance AIC
## <none> 487.36 292.95
## - Bicep 1 531.26 293.35
## - Ankle 1 555.01 295.54
## - Age 1 589.06 298.52
## - Forearm 1 605.07 299.86
## - Abdomen 1 625.19 301.49
## - Hip 1 680.24 305.71
## - Chest 1 952.60 322.55
## - Height 1 1453.11 343.66

Il ne retient plus que Age + Height + Chest + Abdomen + Hip + Ankle + Bicep + Forearm. Il est donc, ce
qui était attendu, plus parcimonieux que celui obtenu avec l’AIC.

2. Avant de construire les modèles ridge et lasso, on transforme le fichier en matrice :

M<-[Link](Weight~.,data=bodyfat)

Dans les modèles ridge et lasso, on cherche des représentations du type :

W eight = β0 + β1 Age + β2 Height + β3 Chest + etc,

avec des protocoles d’estimation des coefficients différents selon la méthode (ridge ou lasso). Modèle ridge :

library(glmnet)

[Link]<-glmnet(M,bodyfat$Weight,alpha=0)
plot([Link],label=T)

27
 15 15 15 15 15 15

2
20
10

5
16
15
13
9
7
1
40
14
6
11
8
12
0

3
Coefficients

−20
−40

0 10 20 30 40 50

L1 Norm
plot([Link],xvar="lambda",label=T,lwd=2)

15 15 15 15 15

2
20
10

5
16
15
13
9
7
10
14
6
11
8
12
0

3
4
Coefficients

−20
−40

2 4 6 8 10

Log Lambda Valida-

tion croisée pour le modèle ridge :

[Link]<-[Link](M,bodyfat$Weight,alpha=0)
plot([Link])

28
 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15

800
Mean−Squared Error

600
400
200
0

2 4 6 8 10

Log(λ)
lam_ridge<-[Link]$[Link]
modelridge<-coef([Link],s=lam_ridge)
modelridge

## 17 x 1 sparse Matrix of class "dgCMatrix"

## s1
## (Intercept) -341.08962404
## (Intercept) .
## Density 22.13293713
## [Link] 0.08478397
## Age -0.15519468
## Height 2.04328737
## Neck 0.59156155
## Chest 0.76521290
## Abdomen 0.33184067
## Waist 0.83816020
## Hip 0.61588135
## Thigh 0.39523946
## Knee 0.21776209
## Ankle 0.97860480
## Bicep 0.59820394
## Forearm 1.21597668
## Wrist 1.51091665

Modèle lasso :

[Link]<-glmnet(M,bodyfat$Weight,alpha=1)
plot([Link],label=T)

29
 0 3 5 5 9 12

5
2.0

16
1.5
Coefficients

7
15
1.0

13
10
14
0.5

8
11
−0.5 0.0

12

0 2 4 6 8 10

L1 Norm
plot([Link],xvar="lambda",label=T,lwd=2)

12 11 10 9 8 5 2

5
2.0

16
1.5
Coefficients

7
15
1.0

13
10
14
0.5

8
11
−0.5 0.0

12

−3 −2 −1 0 1 2 3

Log Lambda

Validation croisée pour le lasso :

[Link]<-[Link](M,bodyfat$Weight,alpha=1)
plot([Link])

30
 12 12 12 11 11 10 10 10 8 8 8 5 5 4 3 2 0

800
Mean−Squared Error

600
400
200
0

−3 −2 −1 0 1 2 3

Log(λ)
lam_lasso<-[Link]$[Link]
modellasso<-coef([Link],s=lam_lasso)
modellasso

## 17 x 1 sparse Matrix of class "dgCMatrix"

## s1
## (Intercept) -328.49891508
## (Intercept) .
## Density .
## [Link] .
## Age -0.13681083
## Height 2.21132932
## Neck 0.01965849
## Chest 1.12412886
## Abdomen 0.50329295
## Waist .
## Hip 0.79829855
## Thigh 0.33252064
## Knee .
## Ankle 0.80412201
## Bicep 0.52931678
## Forearm 1.26351281
## Wrist 1.44077219

3. Nombre de coefficients nuls du lasso :

sum(coef([Link],s=lam_lasso)==0)

## [1] 5

4. On compare 4 modèles par K-fold VC : MCO pas à pas AIC, MCO pas à pas BIC, Ridge, Lasso et MCO.

31
erreur_ridge<-min([Link]$cvm) #erreur CV minimale pour le ridge
lambda_min_ridge<-[Link]$[Link] #lambda correspondant
c(erreur_ridge,lambda_min_ridge)

## [1] 17.764720 2.564556

erreur_lasso<-min([Link]$cvm)
lambda_min_lasso<-[Link]$[Link]
c(erreur_lasso,lambda_min_lasso)

## [1] 14.9314941 0.1687307

library(boot)
[Link](bodyfat,[Link],K=10)$delta[1]

## [1] 14.8921
[Link](bodyfat,[Link],K=10)$delta[1]

## [1] 13.58155
[Link](bodyfat,[Link],K=10)$delta[1]

## [1] 38.07168

5. Pour finir, le modèle elastic net.

library(caret)

## Loading required package: ggplot2

## Loading required package: lattice
##
## Attaching package: 'lattice'
## The following object is masked from 'package:boot':
##
## melanoma
cv10 <- trainControl(method = "cv", number = 10) #Partionnement en vue d'une CV 10-fold
elastic <- train(Weight~.,data = bodyfat,method = "glmnet",trControl = cv10,tuneLength=10)
#Recherche des meilleurs alpha et lambda (chacun sur tuneLength valeurs)
elastic

## glmnet
##
## 50 samples
## 15 predictors
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 46, 46, 46, 46, 45, 46, ...
## Resampling results across tuning parameters:
##
## alpha lambda RMSE Rsquared MAE
## 0.1 0.02737253 3.904359 0.9804478 3.367158
## 0.1 0.06323409 3.904359 0.9804478 3.367158

32
## 0.1 0.14607894 3.868352 0.9809024 3.345703
## 0.1 0.33746130 3.828068 0.9817409 3.331591
## 0.1 0.77957938 3.749877 0.9832398 3.306182
## 0.1 1.80092948 3.854812 0.9832413 3.357714
## 0.1 4.16038067 4.277283 0.9819423 3.555288
## 0.1 9.61101895 5.213984 0.9793252 4.145427
## 0.1 22.20270032 7.141399 0.9758723 6.031587
## 0.2 0.02737253 3.896031 0.9804922 3.360959
## 0.2 0.06323409 3.892998 0.9805245 3.359269
## 0.2 0.14607894 3.846550 0.9811375 3.334057
## 0.2 0.33746130 3.793496 0.9821600 3.307586
## 0.2 0.77957938 3.687161 0.9835208 3.256309
## 0.2 1.80092948 3.829273 0.9830555 3.307118
## 0.2 4.16038067 4.432965 0.9807635 3.571591
## 0.2 9.61101895 5.768060 0.9774802 4.745727
## 0.2 22.20270032 8.959671 0.9733824 7.658441
## 0.3 0.02737253 3.916701 0.9800727 3.377214
## 0.3 0.06323409 3.901564 0.9802669 3.368251
## 0.3 0.14607894 3.840417 0.9811369 3.335778
## 0.3 0.33746130 3.757341 0.9825531 3.293094
## 0.3 0.77957938 3.656746 0.9834380 3.224269
## 0.3 1.80092948 3.855024 0.9825390 3.278259
## 0.3 4.16038067 4.670880 0.9796678 3.721015
## 0.3 9.61101895 6.539520 0.9761740 5.502054
## 0.3 22.20270032 11.115735 0.9678846 9.551251
## 0.4 0.02737253 3.918972 0.9800064 3.380117
## 0.4 0.06323409 3.889458 0.9803657 3.363835
## 0.4 0.14607894 3.830755 0.9813248 3.335581
## 0.4 0.33746130 3.724698 0.9827324 3.269365
## 0.4 0.77957938 3.649001 0.9831183 3.206483
## 0.4 1.80092948 3.924269 0.9817885 3.259287
## 0.4 4.16038067 4.894973 0.9789338 3.909102
## 0.4 9.61101895 7.416488 0.9739181 6.329955
## 0.4 22.20270032 13.138095 0.9587895 11.256306
## 0.5 0.02737253 3.919979 0.9799654 3.381464
## 0.5 0.06323409 3.879747 0.9804565 3.360860
## 0.5 0.14607894 3.828824 0.9814735 3.338629
## 0.5 0.33746130 3.695444 0.9828030 3.244263
## 0.5 0.77957938 3.655324 0.9827403 3.191059
## 0.5 1.80092948 4.000384 0.9810713 3.263151
## 0.5 4.16038067 5.194919 0.9778505 4.216973
## 0.5 9.61101895 8.356053 0.9701378 7.159878
## 0.5 22.20270032 14.960076 0.9555111 12.734570
## 0.6 0.02737253 3.908222 0.9800617 3.378023
## 0.6 0.06323409 3.871811 0.9805749 3.359610
## 0.6 0.14607894 3.820918 0.9816489 3.335964
## 0.6 0.33746130 3.671867 0.9828646 3.224870
## 0.6 0.77957938 3.676463 0.9822845 3.179676
## 0.6 1.80092948 4.053940 0.9805304 3.298228
## 0.6 4.16038067 5.547538 0.9761164 4.561370
## 0.6 9.61101895 9.278227 0.9646831 7.915526
## 0.6 22.20270032 16.850914 0.9521829 14.247626
## 0.7 0.02737253 3.910110 0.9800446 3.378658
## 0.7 0.06323409 3.866980 0.9806364 3.357223

33
## 0.7 0.14607894 3.808232 0.9818270 3.328779
## 0.7 0.33746130 3.656815 0.9828203 3.205680
## 0.7 0.77957938 3.706588 0.9817431 3.181313
## 0.7 1.80092948 4.105995 0.9800012 3.317324
## 0.7 4.16038067 5.906136 0.9737597 4.893765
## 0.7 9.61101895 10.221304 0.9563572 8.656847
## 0.7 22.20270032 18.829817 0.9418374 15.808159
## 0.8 0.02737253 3.906327 0.9800930 3.377049
## 0.8 0.06323409 3.858331 0.9807629 3.354332
## 0.8 0.14607894 3.788598 0.9820260 3.314140
## 0.8 0.33746130 3.653561 0.9826899 3.194090
## 0.8 0.77957938 3.742458 0.9812054 3.177223
## 0.8 1.80092948 4.198364 0.9793658 3.368489
## 0.8 4.16038067 6.245100 0.9711539 5.193617
## 0.8 9.61101895 10.903221 0.9521154 9.196138
## 0.8 22.20270032 20.643973 0.9324841 17.318186
## 0.9 0.02737253 3.910672 0.9800266 3.380447
## 0.9 0.06323409 3.859722 0.9807941 3.356719
## 0.9 0.14607894 3.775378 0.9820944 3.304696
## 0.9 0.33746130 3.658850 0.9824491 3.189699
## 0.9 0.77957938 3.787910 0.9806940 3.184608
## 0.9 1.80092948 4.324150 0.9781847 3.461354
## 0.9 4.16038067 6.545748 0.9682974 5.445116
## 0.9 9.61101895 11.578372 0.9472439 9.785252
## 0.9 22.20270032 22.462527 0.9145548 18.845368
## 1.0 0.02737253 3.908117 0.9800633 3.379075
## 1.0 0.06323409 3.855937 0.9809032 3.355537
## 1.0 0.14607894 3.763085 0.9821783 3.293348
## 1.0 0.33746130 3.665994 0.9822404 3.193003
## 1.0 0.77957938 3.831466 0.9802549 3.186376
## 1.0 1.80092948 4.452385 0.9770545 3.572268
## 1.0 4.16038067 6.783234 0.9669802 5.672091
## 1.0 9.61101895 12.306526 0.9399312 10.384691
## 1.0 22.20270032 24.468825 0.8608662 20.528370
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0.4 and lambda = 0.7795794.

On extrait les meilleures valeurs de alpha et lambda, ainsi que l’erreur par VC 10-fold.

best<-rownames(elastic$bestTune) #Ligne qui donne la plus petite erreur

best_result <- elastic$results[best, ]
c((best_result$RMSE)ˆ2,best_result$alpha,best_result$lambda)

## [1] 13.3152069 0.4000000 0.7795794

6 et 7. Une amélioration possible serait de considérer les variables explicatives retenues à l’issue de la
procédure pas à pas descendante avec l’AIC (par exemple, mais dans une étude complète, il faut aussi
considérer le BIC), et de construire avec ces variables des modèles ridge et lasso. Tout d’abord avec le ridge :

[Link]<-[Link](Weight~Age + Height + Chest + Abdomen + Hip + Thigh + Ankle + Bicep + Forearm + Wris

34
[Link]<-[Link]([Link],bodyfat$Weight,alpha=0)
c(min([Link]$cvm),[Link]$[Link] )

## [1] 14.398742 2.564556

Puis avec le lasso :

[Link]<-[Link]([Link],bodyfat$Weight,alpha=1)
c(min([Link]$cvm),[Link]$[Link] )

## [1] 15.83139259 0.09655388

Enfin, essayons aussi avec elastic net :

[Link] <- train(Weight~Age+Height+Chest+Abdomen+Hip+Thigh+Ankle+Bicep+Forearm+Wrist,data = bodyfat,

[Link]<-rownames([Link]$bestTune)
best_result.AIC <- [Link]$results[[Link], ]
c((best_result.AIC$RMSE)ˆ2,best_result.AIC$alpha,best_result.AIC$lambda)

## [1] 12.5884728 0.1000000 0.7795794

35
EXERCICE 2

Il suffit de prendre l’étude suivie précédemment. On se contente ici de partir du modèle complet, puis de
construire un modèle par procédure pas à pas descendante AIC, et enfin d’en déduire alors les modèles lasso
et ridge.
prostate<-[Link]("~/Dropbox/ABIDJAN/DONNEES/[Link]",header=T)
head(prostate)

## lcavol lweight age lbph svi lcp gleason pgg45 lpsa

## 1 -0.5798185 2.769459 50 -1.386294 0 -1.386294 6 0 -0.4307829
## 2 -0.9942523 3.319626 58 -1.386294 0 -1.386294 6 0 -0.1625189
## 3 -0.5108256 2.691243 74 -1.386294 0 -1.386294 7 20 -0.1625189
## 4 -1.2039728 3.282789 58 -1.386294 0 -1.386294 6 0 -0.1625189
## 5 0.7514161 3.432373 62 -1.386294 0 -1.386294 6 0 0.3715636
## 6 -1.0498221 3.228826 50 -1.386294 0 -1.386294 6 0 0.7654678
[Link]<-glm(lpsa~.,data=prostate,family="gaussian")

[Link]<-step([Link],direction="backward")

## Start: AIC=216.5
## lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason +
## pgg45
##
## Df Deviance AIC
## - gleason 1 43.108 214.61
## - pgg45 1 43.569 215.64
## - lcp 1 43.740 216.02
## <none> 43.058 216.50
## - lbph 1 44.423 217.52
## - age 1 44.857 218.46
## - lweight 1 47.749 224.53
## - svi 1 47.939 224.91
## - lcavol 1 63.258 251.81
##
## Step: AIC=214.61
## lpsa ~ lcavol + lweight + age + lbph + svi + lcp + pgg45
##
## Df Deviance AIC
## - lcp 1 43.776 214.10
## <none> 43.108 214.61
## - pgg45 1 44.306 215.27
## - lbph 1 44.492 215.67
## - age 1 44.865 216.48
## - lweight 1 47.751 222.53
## - svi 1 47.941 222.91
## - lcavol 1 64.427 251.58
##
## Step: AIC=214.1
## lpsa ~ lcavol + lweight + age + lbph + svi + pgg45
##
## Df Deviance AIC
## - pgg45 1 44.437 213.55
## <none> 43.776 214.10

36
## - lbph 1 45.109 215.01
## - age 1 45.264 215.34
## - svi 1 47.953 220.94
## - lweight 1 48.431 221.90
## - lcavol 1 66.531 252.70
##
## Step: AIC=213.55
## lpsa ~ lcavol + lweight + age + lbph + svi
##
## Df Deviance AIC
## <none> 44.437 213.55
## - age 1 45.595 214.05
## - lbph 1 45.945 214.79
## - lweight 1 48.751 220.54
## - svi 1 50.288 223.55
## - lcavol 1 70.379 256.15

La fonction step a donc retenu les variables lcavol + lweight + age + lbph + svi. Examinons les erreurs
K-fold de ces 2 modèles :
[Link](prostate,[Link],K=10)$delta[1]

## [1] 0.5327499
[Link](prostate,[Link],K=10)$delta[1]

## [1] 0.5076091
Le modèle step est le meilleur des 2. On construit maintenant les modèles ridge et lasso avec les variables
retenues par step :
[Link]<-[Link](lpsa~lcavol + lweight + age + lbph + svi,data=prostate)

[Link]<-[Link]([Link],prostate$lpsa,alpha=0)
c(min([Link]$cvm),[Link]$[Link])

## [1] 0.52311635 0.08434274

[Link]<-[Link]([Link],prostate$lpsa,alpha=1)
c(min([Link]$cvm),[Link]$[Link])

## [1] 0.54633546 0.02961435

Le modèle ridge remporte le match. On pourrait également tester des modèles avec interaction.

37
SEANCE DE TP5

EXERCICE 1

[Link]<-[Link]("~/Dropbox/ABIDJAN/DONNEES/[Link]")
head([Link])

## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18

## 1 A11 6 A34 A43 1169 A65 A75 4 A93 A101 4 A121 67 A143 A152 2 A173 1
## 2 A12 48 A32 A43 5951 A61 A73 2 A92 A101 2 A121 22 A143 A152 1 A173 1
## 3 A14 12 A34 A46 2096 A61 A74 2 A93 A101 3 A121 49 A143 A152 1 A172 2
## 4 A11 42 A32 A42 7882 A61 A74 2 A93 A103 4 A122 45 A143 A153 1 A173 2
## 5 A11 24 A33 A40 4870 A61 A73 3 A93 A101 4 A124 53 A143 A153 2 A173 2
## 6 A14 36 A32 A46 9055 A65 A73 2 A93 A101 4 A124 35 A143 A153 1 A172 2
## V19 V20 V21
## 1 A192 A201 1
## 2 A191 A201 2
## 3 A191 A201 1
## 4 A191 A201 1
## 5 A191 A201 2
## 6 A192 A201 1

1. Il faut ramener la variable à expliquer dans {0, 1} et puisqu’on sépare en échantillon test et d’apprentissage,
il faut permuter aléatoirement les données par sécurité. Par convention, 1 si le client est “bon”, 0 sinon.

[Link]$V21<-[Link]$V21 #On ramène la variable à expliquer dans {0,1}

[Link]$V21<-[Link]([Link]$V21)
[Link]<-[Link][sample(nrow([Link])),]
napp<-800
[Link]<-[Link][1:napp,]
[Link]<-[Link][-(1:napp),]

2. Premier modèle, avec toutes les variables.

logitComplet<-glm(V21~.,data=[Link],family=binomial)
summary(logitComplet)

##
## Call:
## glm(formula = V21 ~ ., family = binomial, data = [Link])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0801 -0.6109 0.3490 0.6568 2.0644
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.740e+00 1.294e+00 -2.117 0.034262 *
## V1A12 4.116e-01 2.546e-01 1.617 0.105972
## V1A13 1.205e+00 4.578e-01 2.631 0.008502 **
## V1A14 1.708e+00 2.652e-01 6.438 1.21e-10 ***
## V2 -2.246e-02 1.071e-02 -2.097 0.035993 *

38
## V3A31 2.074e-01 6.780e-01 0.306 0.759736
## V3A32 1.242e+00 5.457e-01 2.276 0.022850 *
## V3A33 1.257e+00 5.973e-01 2.104 0.035370 *
## V3A34 1.748e+00 5.542e-01 3.153 0.001615 **
## V4A41 2.097e+00 4.394e-01 4.773 1.82e-06 ***
## V4A410 1.853e+00 8.132e-01 2.279 0.022661 *
## V4A42 9.787e-01 3.089e-01 3.169 0.001532 **
## V4A43 9.209e-01 2.798e-01 3.292 0.000995 ***
## V4A44 3.482e-01 8.773e-01 0.397 0.691405
## V4A45 4.569e-01 6.634e-01 0.689 0.491019
## V4A46 1.448e-01 4.574e-01 0.317 0.751550
## V4A48 2.159e+00 1.250e+00 1.726 0.084307 .
## V4A49 7.414e-01 3.835e-01 1.933 0.053217 .
## V5 -1.412e-04 5.216e-05 -2.708 0.006772 **
## V6A62 4.266e-01 3.318e-01 1.286 0.198483
## V6A63 4.465e-01 4.415e-01 1.011 0.311853
## V6A64 1.950e+00 7.511e-01 2.596 0.009443 **
## V6A65 8.216e-01 3.039e-01 2.704 0.006854 **
## V7A72 2.028e-01 5.075e-01 0.400 0.689458
## V7A73 3.465e-01 4.827e-01 0.718 0.472816
## V7A74 1.092e+00 5.165e-01 2.114 0.034529 *
## V7A75 5.185e-01 4.865e-01 1.066 0.286517
## V8 -4.450e-01 1.077e-01 -4.130 3.62e-05 ***
## V9A92 7.328e-01 4.466e-01 1.641 0.100830
## V9A93 1.012e+00 4.362e-01 2.320 0.020357 *
## V9A94 8.981e-01 5.226e-01 1.719 0.085705 .
## V10A102 -6.599e-01 5.000e-01 -1.320 0.186867
## V10A103 8.379e-01 4.552e-01 1.841 0.065648 .
## V11 -3.761e-02 1.003e-01 -0.375 0.707679
## V12A122 -4.652e-01 2.919e-01 -1.594 0.111024
## V12A123 -9.836e-02 2.694e-01 -0.365 0.715070
## V12A124 -7.149e-01 4.869e-01 -1.468 0.142021
## V13 2.748e-02 1.150e-02 2.390 0.016843 *
## V14A142 1.030e+00 5.037e-01 2.044 0.040952 *
## V14A143 9.975e-01 2.732e-01 3.651 0.000262 ***
## V15A152 6.301e-01 2.687e-01 2.345 0.019052 *
## V15A153 5.491e-01 5.551e-01 0.989 0.322496
## V16 8.653e-03 2.272e-01 0.038 0.969616
## V17A172 -2.019e-01 8.124e-01 -0.249 0.803690
## V17A173 -3.233e-01 7.769e-01 -0.416 0.677326
## V17A174 -1.623e-01 7.679e-01 -0.211 0.832621
## V18 -2.060e-01 2.866e-01 -0.719 0.472390
## V19A192 3.621e-01 2.353e-01 1.539 0.123760
## V20A202 1.719e+00 6.729e-01 2.554 0.010638 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 968.76 on 799 degrees of freedom
## Residual deviance: 682.38 on 751 degrees of freedom
## AIC: 780.38
##
## Number of Fisher Scoring iterations: 5

39
3. Second modèle, une amélioration de logitComplet par une procédure de descente pas à pas avec le critère
AIC.

logitBAIC<-step(logitComplet,direction="backward")

## Start: AIC=780.38
## V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 +
## V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20
##
## Df Deviance AIC
## - V17 3 682.86 774.86
## - V16 1 682.38 778.38
## - V11 1 682.52 778.52
## - V12 3 686.61 778.61
## - V18 1 682.89 778.89
## - V9 3 688.12 780.12
## <none> 682.38 780.38
## - V19 1 684.77 780.77
## - V7 4 691.48 781.48
## - V15 2 687.89 781.89
## - V10 2 688.09 782.09
## - V2 1 686.77 782.77
## - V13 1 688.33 784.33
## - V5 1 689.85 785.85
## - V20 1 690.43 786.43
## - V6 4 698.61 788.61
## - V3 4 699.00 789.00
## - V14 2 696.04 790.04
## - V8 1 700.33 796.33
## - V4 9 719.59 799.59
## - V1 3 733.28 825.28
##
## Step: AIC=774.86
## V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 +
## V12 + V13 + V14 + V15 + V16 + V18 + V19 + V20
##
## Df Deviance AIC
## - V16 1 682.86 772.86
## - V11 1 683.04 773.04
## - V12 3 687.22 773.22
## - V18 1 683.35 773.35
## - V9 3 688.65 774.65
## <none> 682.86 774.86
## - V19 1 685.60 775.60
## - V7 4 691.65 775.65
## - V15 2 688.29 776.29
## - V10 2 688.51 776.51
## - V2 1 687.86 777.86
## - V13 1 689.41 779.41
## - V5 1 690.11 780.11
## - V20 1 691.09 781.09
## - V6 4 698.93 782.93
## - V3 4 699.44 783.44
## - V14 2 696.47 784.47

40
## - V8 1 700.69 790.69
## - V4 9 720.68 794.68
## - V1 3 733.46 819.46
##
## Step: AIC=772.86
## V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 +
## V12 + V13 + V14 + V15 + V18 + V19 + V20
##
## Df Deviance AIC
## - V11 1 683.04 771.04
## - V12 3 687.24 771.24
## - V18 1 683.35 771.35
## - V9 3 688.66 772.66
## <none> 682.86 772.86
## - V19 1 685.61 773.61
## - V7 4 691.66 773.66
## - V15 2 688.30 774.30
## - V10 2 688.52 774.52
## - V2 1 687.87 775.87
## - V13 1 689.42 777.42
## - V5 1 690.11 778.11
## - V20 1 691.12 779.12
## - V6 4 698.93 780.93
## - V14 2 696.59 782.59
## - V3 4 700.78 782.78
## - V8 1 700.69 788.69
## - V4 9 720.68 792.68
## - V1 3 733.47 817.47
##
## Step: AIC=771.04
## V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V12 +
## V13 + V14 + V15 + V18 + V19 + V20
##
## Df Deviance AIC
## - V12 3 687.39 769.39
## - V18 1 683.54 769.54
## - V9 3 688.80 770.80
## <none> 683.04 771.04
## - V7 4 691.68 771.68
## - V19 1 685.69 771.69
## - V10 2 688.61 772.61
## - V15 2 689.35 773.35
## - V2 1 688.03 774.03
## - V13 1 689.46 775.46
## - V5 1 690.23 776.23
## - V20 1 691.40 777.40
## - V6 4 699.02 779.02
## - V14 2 696.66 780.66
## - V3 4 700.83 780.83
## - V8 1 700.89 786.89
## - V4 9 720.70 790.70
## - V1 3 733.98 815.98
##
## Step: AIC=769.39

41
## V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V13 +
## V14 + V15 + V18 + V19 + V20
##
## Df Deviance AIC
## - V18 1 687.85 767.85
## - V9 3 693.03 769.03
## <none> 687.39 769.39
## - V19 1 689.45 769.45
## - V7 4 696.70 770.70
## - V10 2 693.00 771.00
## - V2 1 692.41 772.41
## - V15 2 695.03 773.03
## - V13 1 693.81 773.81
## - V20 1 695.10 775.10
## - V5 1 695.61 775.61
## - V6 4 702.75 776.75
## - V3 4 705.38 779.38
## - V14 2 701.66 779.66
## - V8 1 706.31 786.31
## - V4 9 725.68 789.68
## - V1 3 741.42 817.42
##
## Step: AIC=767.85
## V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V13 +
## V14 + V15 + V19 + V20
##
## Df Deviance AIC
## - V9 3 693.12 767.12
## <none> 687.85 767.85
## - V19 1 689.92 767.92
## - V7 4 697.11 769.11
## - V10 2 693.39 769.39
## - V2 1 692.84 770.84
## - V15 2 695.46 771.46
## - V13 1 694.09 772.09
## - V20 1 695.53 773.53
## - V5 1 695.89 773.89
## - V6 4 702.98 774.98
## - V3 4 706.25 778.25
## - V14 2 702.45 778.45
## - V8 1 706.41 784.41
## - V4 9 726.03 788.03
## - V1 3 742.12 816.12
##
## Step: AIC=767.12
## V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V10 + V13 + V14 +
## V15 + V19 + V20
##
## Df Deviance AIC
## <none> 693.12 767.12
## - V19 1 695.34 767.34
## - V10 2 699.08 769.08
## - V7 4 703.74 769.74
## - V2 1 698.43 770.43

42
## - V13 1 698.49 770.49
## - V15 2 700.96 770.96
## - V5 1 700.98 772.98
## - V20 1 701.15 773.15
## - V6 4 708.68 774.68
## - V3 4 711.12 777.12
## - V14 2 707.35 777.35
## - V8 1 709.73 781.73
## - V4 9 731.23 787.23
## - V1 3 748.62 816.62

Examinons les caractéristiques de logitBAIC :

summary(logitBAIC)

##
## Call:
## glm(formula = V21 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V10 +
## V13 + V14 + V15 + V19 + V20, family = binomial, data = [Link])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0316 -0.6659 0.3589 0.6630 2.0191
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.546e+00 8.891e-01 -2.863 0.004192 **
## V1A12 4.548e-01 2.493e-01 1.824 0.068091 .
## V1A13 1.290e+00 4.450e-01 2.899 0.003739 **
## V1A14 1.754e+00 2.606e-01 6.731 1.69e-11 ***
## V2 -2.390e-02 1.038e-02 -2.303 0.021283 *
## V3A31 1.745e-01 6.553e-01 0.266 0.789981
## V3A32 1.207e+00 5.300e-01 2.278 0.022742 *
## V3A33 1.203e+00 5.897e-01 2.039 0.041427 *
## V3A34 1.706e+00 5.480e-01 3.113 0.001855 **
## V4A41 2.038e+00 4.287e-01 4.753 2.00e-06 ***
## V4A410 1.887e+00 7.811e-01 2.416 0.015692 *
## V4A42 8.178e-01 2.994e-01 2.731 0.006309 **
## V4A43 9.625e-01 2.746e-01 3.504 0.000458 ***
## V4A44 3.074e-01 8.651e-01 0.355 0.722342
## V4A45 4.529e-01 6.539e-01 0.692 0.488624
## V4A46 3.471e-02 4.499e-01 0.077 0.938499
## V4A48 2.187e+00 1.257e+00 1.740 0.081871 .
## V4A49 7.190e-01 3.740e-01 1.922 0.054564 .
## V5 -1.363e-04 4.902e-05 -2.781 0.005420 **
## V6A62 3.242e-01 3.214e-01 1.009 0.313092
## V6A63 4.658e-01 4.357e-01 1.069 0.284992
## V6A64 1.911e+00 7.280e-01 2.625 0.008659 **
## V6A65 7.650e-01 2.975e-01 2.571 0.010129 *
## V7A72 4.922e-02 4.398e-01 0.112 0.910875
## V7A73 2.621e-01 4.068e-01 0.644 0.519419
## V7A74 1.029e+00 4.542e-01 2.266 0.023456 *
## V7A75 4.026e-01 4.229e-01 0.952 0.341149
## V8 -4.111e-01 1.034e-01 -3.975 7.03e-05 ***

43
## V10A102 -6.172e-01 4.914e-01 -1.256 0.209151
## V10A103 8.459e-01 4.382e-01 1.930 0.053587 .
## V13 2.494e-02 1.096e-02 2.276 0.022846 *
## V14A142 1.122e+00 4.971e-01 2.258 0.023934 *
## V14A143 9.826e-01 2.660e-01 3.694 0.000221 ***
## V15A152 6.441e-01 2.514e-01 2.562 0.010411 *
## V15A153 1.493e-01 4.056e-01 0.368 0.712808
## V19A192 3.266e-01 2.202e-01 1.483 0.138108
## V20A202 1.651e+00 6.482e-01 2.548 0.010841 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 968.76 on 799 degrees of freedom
## Residual deviance: 693.12 on 763 degrees of freedom
## AIC: 767.12
##
## Number of Fisher Scoring iterations: 5

Puis une procédure pas à pas ascendante avec le critère AIC :

model0<-glm(V21~1,data=[Link],family="binomial")
logitFAIC<-step(model0,scope=formula(logitComplet),direction="forward")

## Start: AIC=970.76
## V21 ~ 1
##
## Df Deviance AIC
## + V1 3 868.14 876.14
## + V3 4 919.60 929.60
## + V2 1 935.58 939.58
## + V6 4 936.60 946.60
## + V12 3 949.07 957.07
## + V4 9 937.31 957.31
## + V15 2 951.32 957.32
## + V5 1 954.74 958.74
## + V13 1 955.42 959.42
## + V14 2 953.83 959.83
## + V7 4 953.99 963.99
## + V8 1 960.32 964.32
## + V20 1 962.03 966.03
## + V16 1 965.02 969.02
## + V10 2 964.33 970.33
## + V9 3 962.46 970.46
## <none> 968.76 970.76
## + V19 1 966.78 970.78
## + V18 1 968.76 972.76
## + V11 1 968.76 972.76
## + V17 3 967.09 975.09
##
## Step: AIC=876.14
## V21 ~ V1
##

44
## Df Deviance AIC
## + V2 1 840.16 850.16
## + V3 4 836.34 852.34
## + V12 3 851.06 865.06
## + V4 9 839.49 865.49
## + V5 1 855.91 865.91
## + V14 2 854.70 866.70
## + V8 1 858.79 868.79
## + V6 4 853.42 869.42
## + V20 1 859.49 869.49
## + V13 1 859.59 869.59
## + V15 2 857.84 869.84
## + V10 2 859.54 871.54
## + V7 4 857.39 873.39
## <none> 868.14 876.14
## + V16 1 866.45 876.45
## + V19 1 867.62 877.62
## + V11 1 867.99 877.99
## + V9 3 864.06 878.06
## + V18 1 868.14 878.14
## + V17 3 865.49 879.49
##
## Step: AIC=850.16
## V21 ~ V1 + V2
##
## Df Deviance AIC
## + V3 4 814.56 832.56
## + V4 9 805.29 833.29
## + V14 2 827.11 841.11
## + V7 4 824.11 842.11
## + V6 4 824.82 842.82
## + V13 1 831.16 843.16
## + V8 1 832.57 844.57
## + V15 2 832.43 846.43
## + V20 1 834.60 846.60
## + V10 2 833.01 847.01
## + V12 3 833.69 849.69
## + V19 1 837.87 849.87
## <none> 840.16 850.16
## + V9 3 834.67 850.67
## + V16 1 838.81 850.81
## + V11 1 839.77 851.77
## + V5 1 840.03 852.03
## + V18 1 840.15 852.15
## + V17 3 839.78 855.78
##
## Step: AIC=832.56
## V21 ~ V1 + V2 + V3
##
## Df Deviance AIC
## + V4 9 778.64 814.64
## + V6 4 799.27 825.27
## + V8 1 805.45 825.45
## + V7 4 801.17 827.17

45
## + V13 1 807.41 827.41
## + V14 2 805.47 827.47
## + V10 2 806.91 828.91
## + V20 1 809.07 829.07
## + V15 2 809.14 831.14
## + V19 1 812.04 832.04
## <none> 814.56 832.56
## + V12 3 809.11 833.11
## + V9 3 809.45 833.45
## + V18 1 814.37 834.37
## + V11 1 814.44 834.44
## + V5 1 814.54 834.54
## + V16 1 814.55 834.55
## + V17 3 814.16 838.16
##
## Step: AIC=814.64
## V21 ~ V1 + V2 + V3 + V4
##
## Df Deviance AIC
## + V6 4 762.54 806.54
## + V8 1 770.67 808.67
## + V20 1 771.40 809.40
## + V14 2 769.56 809.56
## + V7 4 765.95 809.95
## + V13 1 772.18 810.18
## + V15 2 771.66 811.66
## + V10 2 772.46 812.46
## + V12 3 772.12 814.12
## <none> 778.64 814.64
## + V19 1 776.94 814.94
## + V5 1 777.92 815.92
## + V9 3 774.15 816.15
## + V18 1 778.43 816.43
## + V16 1 778.62 816.62
## + V11 1 778.64 816.64
## + V17 3 777.74 819.74
##
## Step: AIC=806.54
## V21 ~ V1 + V2 + V3 + V4 + V6
##
## Df Deviance AIC
## + V8 1 753.97 799.97
## + V20 1 755.54 801.54
## + V14 2 753.73 801.73
## + V7 4 750.60 802.60
## + V10 2 754.82 802.82
## + V13 1 756.84 802.84
## + V15 2 754.98 802.98
## + V12 3 755.47 805.47
## <none> 762.54 806.54
## + V19 1 761.41 807.41
## + V9 3 757.75 807.75
## + V5 1 761.76 807.76
## + V18 1 762.44 808.44

46
## + V11 1 762.48 808.48
## + V16 1 762.52 808.52
## + V17 3 761.50 811.50
##
## Step: AIC=799.97
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8
##
## Df Deviance AIC
## + V14 2 744.10 794.10
## + V13 1 747.34 795.34
## + V15 2 745.47 795.47
## + V7 4 741.73 795.73
## + V20 1 747.96 795.96
## + V5 1 748.43 796.43
## + V10 2 746.57 796.57
## + V9 3 747.01 799.01
## + V12 3 747.60 799.60
## <none> 753.97 799.97
## + V19 1 752.82 800.82
## + V18 1 753.93 801.93
## + V11 1 753.95 801.95
## + V16 1 753.97 801.97
## + V17 3 753.10 805.10
##
## Step: AIC=794.1
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14
##
## Df Deviance AIC
## + V13 1 736.27 788.27
## + V7 4 730.44 788.44
## + V10 2 734.61 788.61
## + V20 1 737.39 789.39
## + V15 2 736.03 790.03
## + V5 1 738.23 790.23
## + V9 3 736.57 792.57
## <none> 744.10 794.10
## + V12 3 738.82 794.82
## + V19 1 743.08 795.08
## + V18 1 743.91 795.91
## + V16 1 744.04 796.04
## + V11 1 744.09 796.09
## + V17 3 742.99 798.99
##
## Step: AIC=788.27
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13
##
## Df Deviance AIC
## + V20 1 728.73 782.73
## + V10 2 726.86 782.86
## + V5 1 729.35 783.35
## + V7 4 723.80 783.80
## + V15 2 728.25 784.25
## + V12 3 727.73 785.73
## + V9 3 728.85 786.85

47
## <none> 736.27 788.27
## + V11 1 735.67 789.67
## + V19 1 736.01 790.01
## + V18 1 736.26 790.26
## + V16 1 736.26 790.26
## + V17 3 735.05 793.05
##
## Step: AIC=782.73
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13 + V20
##
## Df Deviance AIC
## + V5 1 721.27 777.27
## + V10 2 719.88 777.88
## + V7 4 716.55 778.55
## + V15 2 720.73 778.73
## + V12 3 720.54 780.54
## + V9 3 721.91 781.91
## <none> 728.73 782.73
## + V19 1 728.19 784.19
## + V11 1 728.34 784.34
## + V16 1 728.68 784.68
## + V18 1 728.73 784.73
## + V17 3 727.89 787.89
##
## Step: AIC=777.27
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13 + V20 + V5
##
## Df Deviance AIC
## + V7 4 708.70 772.70
## + V15 2 713.62 773.62
## + V10 2 713.88 773.88
## + V9 3 713.72 775.72
## + V19 1 719.14 777.14
## + V12 3 715.20 777.20
## <none> 721.27 777.27
## + V11 1 720.78 778.78
## + V16 1 721.24 779.24
## + V18 1 721.27 779.27
## + V17 3 720.89 782.89
##
## Step: AIC=772.7
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13 + V20 + V5 + V7
##
## Df Deviance AIC
## + V15 2 701.52 769.52
## + V10 2 702.70 770.70
## + V9 3 702.64 772.64
## <none> 708.70 772.70
## + V19 1 706.75 772.75
## + V11 1 707.76 773.76
## + V12 3 703.79 773.79
## + V18 1 708.68 774.68
## + V16 1 708.70 774.70
## + V17 3 708.10 778.10

48
##
## Step: AIC=769.52
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13 + V20 + V5 + V7 +
## V15
##
## Df Deviance AIC
## + V10 2 695.34 767.34
## + V19 1 699.08 769.08
## <none> 701.52 769.52
## + V9 3 695.65 769.65
## + V18 1 701.46 771.46
## + V11 1 701.49 771.49
## + V16 1 701.50 771.50
## + V12 3 697.81 771.81
## + V17 3 700.94 774.94
##
## Step: AIC=767.34
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13 + V20 + V5 + V7 +
## V15 + V10
##
## Df Deviance AIC
## + V19 1 693.12 767.12
## <none> 695.34 767.34
## + V9 3 689.92 767.92
## + V18 1 695.25 769.25
## + V11 1 695.27 769.27
## + V16 1 695.31 769.31
## + V12 3 691.71 769.71
## + V17 3 694.58 772.58
##
## Step: AIC=767.12
## V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13 + V20 + V5 + V7 +
## V15 + V10 + V19
##
## Df Deviance AIC
## <none> 693.12 767.12
## + V9 3 687.85 767.85
## + V12 3 688.88 768.88
## + V11 1 692.99 768.99
## + V18 1 693.03 769.03
## + V16 1 693.10 769.10
## + V17 3 692.41 772.41
summary(logitFAIC)

##
## Call:
## glm(formula = V21 ~ V1 + V2 + V3 + V4 + V6 + V8 + V14 + V13 +
## V20 + V5 + V7 + V15 + V10 + V19, family = "binomial", data = [Link])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0316 -0.6659 0.3589 0.6630 2.0191
##
## Coefficients:

49
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.546e+00 8.891e-01 -2.863 0.004192 **
## V1A12 4.548e-01 2.493e-01 1.824 0.068091 .
## V1A13 1.290e+00 4.450e-01 2.899 0.003739 **
## V1A14 1.754e+00 2.606e-01 6.731 1.69e-11 ***
## V2 -2.390e-02 1.038e-02 -2.303 0.021283 *
## V3A31 1.745e-01 6.553e-01 0.266 0.789981
## V3A32 1.207e+00 5.300e-01 2.278 0.022742 *
## V3A33 1.203e+00 5.897e-01 2.039 0.041427 *
## V3A34 1.706e+00 5.480e-01 3.113 0.001855 **
## V4A41 2.038e+00 4.287e-01 4.753 2.00e-06 ***
## V4A410 1.887e+00 7.811e-01 2.416 0.015692 *
## V4A42 8.178e-01 2.994e-01 2.731 0.006309 **
## V4A43 9.625e-01 2.746e-01 3.504 0.000458 ***
## V4A44 3.074e-01 8.651e-01 0.355 0.722342
## V4A45 4.529e-01 6.539e-01 0.692 0.488624
## V4A46 3.471e-02 4.499e-01 0.077 0.938499
## V4A48 2.187e+00 1.257e+00 1.740 0.081871 .
## V4A49 7.190e-01 3.740e-01 1.922 0.054564 .
## V6A62 3.242e-01 3.214e-01 1.009 0.313092
## V6A63 4.658e-01 4.357e-01 1.069 0.284992
## V6A64 1.911e+00 7.280e-01 2.625 0.008659 **
## V6A65 7.650e-01 2.975e-01 2.571 0.010129 *
## V8 -4.111e-01 1.034e-01 -3.975 7.03e-05 ***
## V14A142 1.122e+00 4.971e-01 2.258 0.023934 *
## V14A143 9.826e-01 2.660e-01 3.694 0.000221 ***
## V13 2.494e-02 1.096e-02 2.276 0.022846 *
## V20A202 1.651e+00 6.482e-01 2.548 0.010841 *
## V5 -1.363e-04 4.902e-05 -2.781 0.005420 **
## V7A72 4.922e-02 4.398e-01 0.112 0.910875
## V7A73 2.621e-01 4.068e-01 0.644 0.519419
## V7A74 1.029e+00 4.542e-01 2.266 0.023456 *
## V7A75 4.026e-01 4.229e-01 0.952 0.341149
## V15A152 6.441e-01 2.514e-01 2.562 0.010411 *
## V15A153 1.493e-01 4.056e-01 0.368 0.712808
## V10A102 -6.172e-01 4.914e-01 -1.256 0.209151
## V10A103 8.459e-01 4.382e-01 1.930 0.053587 .
## V19A192 3.266e-01 2.202e-01 1.483 0.138108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 968.76 on 799 degrees of freedom
## Residual deviance: 693.12 on 763 degrees of freedom
## AIC: 767.12
##
## Number of Fisher Scoring iterations: 5

4. On calcule les erreurs VC K-fold pour chaque modèle :

library(boot)
[Link]([Link],logitComplet,K=10)$delta[1]

50
## [1] 0.1675838
[Link]([Link],logitBAIC,K=10)$delta[1]

## [1] 0.1593015
[Link]([Link],logitFAIC,K=10)$delta[1]

## [1] 0.1601075

5. Calculons, pour les 4 modèles, les prévisions de chaque individu de l’échantillon test.

prevComplet<-predict(logitComplet,newdata=[Link],type="response")
prevBAIC<-predict(logitBAIC,newdata=[Link],type="response")
prevFAIC<-predict(logitFAIC,newdata=[Link],type="response")

Puis on calcule les erreurs pour chaque modèle :

mean(round(prevComplet)!=[Link]$V21)

## [1] 0.285
mean(round(prevBAIC)!=[Link]$V21)

## [1] 0.29
mean(round(prevFAIC)!=[Link]$V21)

## [1] 0.29

6. On peut aussi dresser les tables de confusion pour chaque modèle :

table(prev=round(prevComplet),obs=[Link]$V21)

## obs
## prev 0 1
## 0 29 21
## 1 36 114
table(prev=round(prevBAIC),obs=[Link]$V21)

## obs
## prev 0 1
## 0 30 23
## 1 35 112
table(prev=round(prevFAIC),obs=[Link]$V21)

## obs
## prev 0 1
## 0 30 23
## 1 35 112

7. On trace les courbes ROC des modèles :

51
library(pROC)

## Type 'citation("pROC")' for a citation.

##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
roc([Link]$V21,prevComplet,plot=T)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases
##
## Call:
## [Link](response = [Link]$V21, predictor = prevComplet, plot = T)
##
## Data: prevComplet in 65 controls ([Link]$V21 0) < 135 cases ([Link]$V21 1).
## Area under the curve: 0.751
roc([Link]$V21,prevBAIC,plot=T,col="red",add=T)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases
##
## Call:
## [Link](response = [Link]$V21, predictor = prevBAIC, plot = T, col = "red", add = T)
##
## Data: prevBAIC in 65 controls ([Link]$V21 0) < 135 cases ([Link]$V21 1).
## Area under the curve: 0.7486
roc([Link]$V21,prevFAIC,plot=T,col="green",add=T)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

52
 1.0
0.8
0.6
Sensitivity
0.4
0.2
0.0

1.0 0.5 0.0

Specificity
##
## Call:
## [Link](response = [Link]$V21, predictor = prevFAIC, plot = T, col = "green", add = T)
##
## Data: prevFAIC in 65 controls ([Link]$V21 0) < 135 cases ([Link]$V21 1).
## Area under the curve: 0.7486

Examinons les AUC (déjà calculées au dessus d’ailleurs !) :

auc([Link]$V21,prevComplet)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases
## Area under the curve: 0.751
auc([Link]$V21,prevBAIC)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases
## Area under the curve: 0.7486
auc([Link]$V21,prevFAIC)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases
## Area under the curve: 0.7486

8. Enfin, calculons les odds ratio pour le modèle logitFAIC.

53
library(questionr)
[Link](logitFAIC)

## Waiting for profiling to be done...

## OR 2.5 % 97.5 % p
## (Intercept) 0.078404 0.013339 0.4393 0.0041916 **
## V1A12 1.575862 0.967932 2.5755 0.0680908 .
## V1A13 3.633177 1.569150 9.1078 0.0037393 **
## V1A14 5.775612 3.494271 9.7215 1.691e-11 ***
## V2 0.976381 0.956580 0.9964 0.0212828 *
## V3A31 1.190689 0.330670 4.3737 0.7899813
## V3A32 3.343843 1.209068 9.7918 0.0227417 *
## V3A33 3.328506 1.071711 10.9424 0.0414266 *
## V3A34 5.505040 1.920173 16.6692 0.0018546 **
## V4A41 7.672040 3.416289 18.4834 2.002e-06 ***
## V4A410 6.601172 1.490995 32.9907 0.0156916 *
## V4A42 2.265608 1.266517 4.1042 0.0063095 **
## V4A43 2.618158 1.534421 4.5107 0.0004576 ***
## V4A44 1.359867 0.254215 8.1891 0.7223423
## V4A45 1.572789 0.452107 6.0786 0.4886244
## V4A46 1.035319 0.429515 2.5212 0.9384994
## V4A48 8.904835 1.003639 206.0353 0.0818712 .
## V4A49 2.052395 0.993088 4.3175 0.0545640 .
## V6A62 1.382961 0.742757 2.6276 0.3130924
## V6A63 1.593299 0.703386 3.9332 0.2849916
## V6A64 6.760140 1.886356 35.1568 0.0086590 **
## V6A65 2.148932 1.214652 3.9107 0.0101292 *
## V8 0.662926 0.539469 0.8096 7.029e-05 ***
## V14A142 3.072333 1.181542 8.3727 0.0239340 *
## V14A143 2.671286 1.585928 4.5084 0.0002206 ***
## V13 1.025254 1.003799 1.0479 0.0228461 *
## V20A202 5.214378 1.610532 21.3455 0.0108412 *
## V5 0.999864 0.999766 1.0000 0.0054204 **
## V7A72 1.050456 0.442375 2.4927 0.9108752
## V7A73 1.299613 0.583274 2.8893 0.5194186
## V7A74 2.798556 1.152030 6.8711 0.0234564 *
## V7A75 1.495712 0.649423 3.4243 0.3411491
## V15A152 1.904211 1.163260 3.1220 0.0104108 *
## V15A153 1.161033 0.524446 2.5800 0.7128080
## V10A102 0.539473 0.205561 1.4290 0.2091512
## V10A103 2.329989 1.020510 5.7581 0.0535874 .
## V19A192 1.386220 0.902401 2.1427 0.1381076
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On voit ainsi quelles modalités de quelles variables explicatives ont le plus d’influence sur le succès (qui est
ici le fait que le client ne soit pas “bon” pour les critères de la société).

54
EXERCICE 2

panne<-[Link]("~/Dropbox/ABIDJAN/DONNEES/[Link]",header=T)
head(panne)

## etat age marque

## 1 0 4 A
## 2 0 2 C
## 3 0 3 C
## 4 0 9 B
## 5 0 7 B
## 6 0 6 A
marque est un caractère, on le transforme donc en facteur qui, lui, possède des niveaux.
panne$marque<-factor(panne$marque)
head(panne)

## etat age marque

## 1 0 4 A
## 2 0 2 C
## 3 0 3 C
## 4 0 9 B
## 5 0 7 B
## 6 0 6 A
On cherche tout d’abord à estimer le modèle logistique défini par

logit p(age, marque) = β0 + β1 1A (marque) + β2 1B (marque) + β3 1C (marque) + β4 age,

avec la contrainte d’identifiabilité β1 = 0.

model3<-glm(etat~.,data=panne,family=binomial)
model3

##
## Call: glm(formula = etat ~ ., family = binomial, data = panne)
##
## Coefficients:
## (Intercept) age marqueB marqueC
## 0.47808 0.01388 -0.41941 -1.45608
##
## Degrees of Freedom: 32 Total (i.e. Null); 29 Residual
## Null Deviance: 45.72
## Residual Deviance: 43.5 AIC: 51.5
summary(model3)

##
## Call:
## glm(formula = etat ~ ., family = binomial, data = panne)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4232 -1.2263 0.9082 1.1062 1.5982
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)

55
## (Intercept) 0.47808 0.83301 0.574 0.566
## age 0.01388 0.09398 0.148 0.883
## marqueB -0.41941 0.81428 -0.515 0.607
## marqueC -1.45608 1.05358 -1.382 0.167
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 45.717 on 32 degrees of freedom
## Residual deviance: 43.502 on 29 degrees of freedom
## AIC: 51.502
##
## Number of Fisher Scoring iterations: 4

Au niveau 5%, on ne rejette donc pas les hypothèses de nullité des coefficients, ce qui signifie que compte tenu
du modèle proposé, la meilleure configuration semble être celle où logit pβ (x) = β0 ! En particulier, l’âge et
la marque n’interviendraient pas dans les pannes, ce qui est manifestement en contradiction avec les données
! Le modèle 3 n’ajuste pas bien les données et ne semble donc pas pertinent. Pourquoi ? On commence par
étudier dans les questions qui suivent le modèle selon age uniquement.

2. On considère tout d’abord la seule variable explicative age :

model5<-glm(etat~age,data=panne,family=binomial)
model5

##
## Call: glm(formula = etat ~ age, family = binomial, data = panne)
##
## Coefficients:
## (Intercept) age
## -0.10748 0.03141
##
## Degrees of Freedom: 32 Total (i.e. Null); 31 Residual
## Null Deviance: 45.72
## Residual Deviance: 45.6 AIC: 49.6
Le modèle estimé est donc, pour la variable age :
1
pb(age) = .
1 + e0.11−0.03∗age

Ici, x = age.
c5<-model5$coefficients #Extraction des coefficients de model5
f5=function(x) 1/(1+exp(-c5[1]-c5[2]*x))
curve(expr=f5, main="Probabilité estimée", from=1, to= 16)

56
 Probabilité estimée
0.60
0.56
f5(x)

0.52
0.48

5 10 15

x
Avant de dessiner le graphe, rentrons les données dans un tableau :
tab<-table(panne$etat,panne$age)
tab

##
## 1 2 3 4 5 6 7 8 9 10 13 14 16
## 0 0 2 2 3 2 3 2 0 2 0 0 0 0
## 1 5 2 2 0 2 0 0 1 1 1 1 1 1
La table tab compte les 0 et les 1 (colonne) par âge (ligne). Graphiquement, cela donne en terme de fréquence
:
freq<-tab[2,]/apply(tab,2,sum)
plot([Link](colnames(tab)),freq, main="Fréquence des pannes en fonction de l'âge")

57
 Fréquence des pannes en fonction de l'âge
1.0
0.8
0.6
freq

0.4
0.2
0.0

5 10 15

[Link](colnames(tab))
En comparaison avec le graphique plus haut, ce n’est pas du tout convaincant.

3. On introduit maintenant la variable explicative “carré de l’âge”.

a2<-panne$ageˆ2
model6<-glm(etat~age+a2,data=panne,family=binomial)
summary(model6)

##
## Call:
## glm(formula = etat ~ age + a2, family = binomial, data = panne)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.54043 -0.74739 0.00033 0.64877 1.88091
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.18501 1.73860 2.407 0.01608 *
## age -2.03343 0.77401 -2.627 0.00861 **
## a2 0.17601 0.07044 2.499 0.01247 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 45.717 on 32 degrees of freedom
## Residual deviance: 31.279 on 30 degrees of freedom
## AIC: 37.279
##

58
## Number of Fisher Scoring iterations: 6
Le modèle devient
1
pb(age) = .
1+ e−4.18+2.03∗age−0.17∗age2

Pour comparer les 2 modèles, on peut utiliser des critères d’information du type AIC (Akaike Information
Criterion) ou BIC (Bayesian Information Criterion).

Comparons model5 et model6 en terme d’AIC et de BIC :

c(AIC(model5),AIC(model6))

## [1] 49.5978 37.2792

c(BIC(model5),BIC(model6))

## [1] 52.59081 41.76872

Model6 est donc choisi par ces 2 critères. On peut aussi comparer ces 2 modèles emboîtés avec un test de
rapport de vraisemblance. Il donne la même conclusion au niveau 5%.
library(car)

## Loading required package: carData

##
## Attaching package: 'car'
## The following object is masked from 'package:boot':
##
## logit
anova(model5,model6,test="LRT")

## Analysis of Deviance Table

##
## Model 1: etat ~ age
## Model 2: etat ~ age + a2
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 31 45.598
## 2 30 31.279 1 14.319 0.0001543 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Représentons maintenant la probabilité estimée de tomber en panne en fonction de l’âge pour le modèle 6.
c6<-model6$coefficients
f6=function(x) 1/(1+exp(-c6[1]-c6[2]*x-c6[3]*xˆ2))
curve(expr=f6, main="Probabilité estimée", from=1, to= 16)

59
 Probabilité estimée
1.0
0.8
0.6
f6(x)

0.4
0.2

5 10 15

x
Beaucoup plus proche du graphe des données. Et aussi plus cohérent avec ce qu’on attend, i.e. un graphe en
forme de baignoire !

On essaie avec le cube de l’âge en plus.

a3<-panne$ageˆ3
model7<-glm(etat~age+a2+a3,data=panne,family=binomial)
summary(model7)

##
## Call:
## glm(formula = etat ~ age + a2 + a3, family = binomial, data = panne)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.53342 -0.72402 0.00937 0.71691 1.88339
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.731972 2.799311 1.690 0.091 .
## age -2.491277 1.920481 -1.297 0.195
## a2 0.276863 0.381626 0.725 0.468
## a3 -0.006278 0.022743 -0.276 0.783
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 45.717 on 32 degrees of freedom
## Residual deviance: 31.219 on 29 degrees of freedom

60
## AIC: 39.219
##
## Number of Fisher Scoring iterations: 8
Noter que les hypothèses de nullité des coefficients de model7 ne sont pas rejetées.

En comparant les AIC et BIC, on voit qu’il semble inutile d’introduire age3 .
c(AIC(model6),AIC(model7))

## [1] 37.27920 39.21892

c(BIC(model6),BIC(model7))

## [1] 41.76872 45.20495

On peut également faire des tests de vraisemblance entre les 2 modèles emboîtés model6 et model7 et parvenir
à la même conclusion au niveau 5% :
library(car)
anova(model6,model7,test="LRT")

## Analysis of Deviance Table

##
## Model 1: etat ~ age + a2
## Model 2: etat ~ age + a2 + a3
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 30 31.279
## 2 29 31.219 1 0.060284 0.806
La nullité du coefficient supplémentaire a3 de model7 n’est donc pas refusée au niveau 5%.

Voyons sur le graphique de la probabilité estimée (avec x la variable age), le changement n’est vraiment pas
notable !
c7<-model7$coefficients
x<-seq(1,16,by=0.001)
xbet<-c7[1]+c7[2]*x+c7[3]*xˆ2+c7[4]*xˆ3
y<-exp(xbet)/(1+exp(xbet))
plot(x,y,type="l")

61
 1.0
0.8
0.6
y

0.4
0.2

5 10 15

x
4. On introduit maintenant le carré de l’âge dans le premier modèle.
model8<-glm(etat~age+a2+marque,data=panne,family=binomial)
summary(model8)

##
## Call:
## glm(formula = etat ~ age + a2 + marque, family = binomial, data = panne)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.63397 -0.71853 0.00029 0.55197 1.81772
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.67773 2.04692 2.285 0.0223 *
## age -2.07273 0.82361 -2.517 0.0118 *
## a2 0.17854 0.07485 2.385 0.0171 *
## marqueB -0.21690 1.07759 -0.201 0.8405
## marqueC -1.09297 1.26145 -0.866 0.3863
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 45.717 on 32 degrees of freedom
## Residual deviance: 30.382 on 28 degrees of freedom
## AIC: 40.382
##
## Number of Fisher Scoring iterations: 7
On compare maintenant model3 et model8, model3 étant un sous-modèle de model8 :

62
anova(model3,model8,test="LRT")

## Analysis of Deviance Table

##
## Model 1: etat ~ age + marque
## Model 2: etat ~ age + a2 + marque
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 29 43.502
## 2 28 30.382 1 13.12 0.0002922 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
On est donc amené à rejetter au niveau 5% l’hypothèse de nullité des coefficients supplémentaires de model8.
Avec les AIC et BIC, même conclusion :
c(AIC(model3),AIC(model8))

## [1] 51.50218 40.38216

c(BIC(model3),BIC(model8))

## [1] 57.48821 47.86470

Avec model8, on a obtenu les coefficients suivants :
model8$coefficient

## (Intercept) age a2 marqueB marqueC

## 4.6777313 -2.0727299 0.1785412 -0.2169039 -1.0929691

63
EXERCICE 3

On veut expliquer les chd (cardiac diseases) par des variables ou facteurs explicatifs, comme le tabac,
l’adiposité, l’alcool, l’âge, l’historique familial. . .

data(SAheart,package="bestglm")
head(SAheart)

## sbp tobacco ldl adiposity famhist typea obesity alcohol age chd
## 1 160 12.00 5.73 23.11 Present 49 25.30 97.20 52 1
## 2 144 0.01 4.41 28.61 Absent 55 28.87 2.06 63 1
## 3 118 0.08 3.48 32.28 Present 52 29.14 3.81 46 0
## 4 170 7.50 6.41 38.03 Present 51 31.99 24.26 58 1
## 5 134 13.60 3.50 27.78 Present 60 25.99 57.34 49 1
## 6 132 6.20 6.47 36.21 Present 62 30.77 14.14 45 0

1. On construit tout d’abord le modèle avec toutes les variables explicatives.

modelbase<-glm(chd~.,data=SAheart,family=binomial)
summary(modelbase)

##
## Call:
## glm(formula = chd ~ ., family = binomial, data = SAheart)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7781 -0.8213 -0.4387 0.8889 2.5435
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.1507209 1.3082600 -4.701 2.58e-06 ***
## sbp 0.0065040 0.0057304 1.135 0.256374
## tobacco 0.0793764 0.0266028 2.984 0.002847 **
## ldl 0.1739239 0.0596617 2.915 0.003555 **
## adiposity 0.0185866 0.0292894 0.635 0.525700
## famhistPresent 0.9253704 0.2278940 4.061 4.90e-05 ***
## typea 0.0395950 0.0123202 3.214 0.001310 **
## obesity -0.0629099 0.0442477 -1.422 0.155095
## alcohol 0.0001217 0.0044832 0.027 0.978350
## age 0.0452253 0.0121298 3.728 0.000193 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 596.11 on 461 degrees of freedom
## Residual deviance: 472.14 on 452 degrees of freedom
## AIC: 492.14
##
## Number of Fisher Scoring iterations: 5

64
Lançons une procédure de recherche exhaustive pour les critères AIC et BIC (on peut le faire car il n’y pas
trop de variables explicatives) :

library(bestglm)

[Link]<-bestglm(SAheart,family=binomial,IC="BIC")

## Morgan-Tatar search since family is non-gaussian.

[Link]<-bestglm(SAheart,family=binomial,IC="AIC")

## Morgan-Tatar search since family is non-gaussian.

[Link]$BestModel

##
## Call: glm(formula = y ~ ., family = family, data = Xi, weights = weights)
##
## Coefficients:
## (Intercept) tobacco ldl famhistPresent typea
## -6.44644 0.08038 0.16199 0.90818 0.03712
## age
## 0.05046
##
## Degrees of Freedom: 461 Total (i.e. Null); 456 Residual
## Null Deviance: 596.1
## Residual Deviance: 475.7 AIC: 487.7
[Link]$BestModel

##
## Call: glm(formula = y ~ ., family = family, data = Xi, weights = weights)
##
## Coefficients:
## (Intercept) tobacco ldl famhistPresent typea
## -6.44644 0.08038 0.16199 0.90818 0.03712
## age
## 0.05046
##
## Degrees of Freedom: 461 Total (i.e. Null); 456 Residual
## Null Deviance: 596.1
## Residual Deviance: 475.7 AIC: 487.7

Ce sont les mêmes modèles. Ajoutons l’interaction age : famhist dans un modèle réduit tenant compte de
tobacco, ldl, famhist, typea et age :

model.inter1<-glm(chd~tobacco+ldl+famhist+typea+age:famhist,data=SAheart,family=binomial)
model.inter1

##
## Call: glm(formula = chd ~ tobacco + ldl + famhist + typea + age:famhist,
## family = binomial, data = SAheart)
##
## Coefficients:
## (Intercept) tobacco ldl famhistPresent
## -6.04259 0.08373 0.16416 -0.22224

65
## typea famhistAbsent:age famhistPresent:age
## 0.03777 0.04041 0.06429
##
## Degrees of Freedom: 461 Total (i.e. Null); 455 Residual
## Null Deviance: 596.1
## Residual Deviance: 474.1 AIC: 488.1
AIC(model.inter1)

## [1] 488.1314

Faisons de même avec le rajout de l’interaction ldl : famhist :

model.inter2<-glm(chd~tobacco+ldl+famhist+typea+age+ldl:famhist,data=SAheart,family=binomial)
model.inter2

##
## Call: glm(formula = chd ~ tobacco + ldl + famhist + typea + age + ldl:famhist,
## family = binomial, data = SAheart)
##
## Coefficients:
## (Intercept) tobacco ldl famhistPresent
## -5.79224 0.08496 0.01758 -0.77068
## typea age ldl:famhistPresent
## 0.03690 0.05140 0.33334
##
## Degrees of Freedom: 461 Total (i.e. Null); 455 Residual
## Null Deviance: 596.1
## Residual Deviance: 466.9 AIC: 480.9
AIC(model.inter2)

## [1] 480.8952

A l’issue de cette étape de l’étude, on retient le modèle model.inter2 qui a le plus faible AIC. Nous essayons
donc de l’améliorer avec une procédure ridge ou lasso.

library(boot)
library(glmnet)
[Link](SAheart,modelbase,K=10)$delta[1]

## [1] 0.1804277
[Link](SAheart,model.inter2,K=10)$delta[1]

## [1] 0.1747501
[Link]<-[Link](chd~tobacco + ldl + famhist + typea + age + ldl:famhist,data=SAheart)[,-1]

[Link]<-[Link]([Link],SAheart$chd,alpha=0)
min([Link]$cvm)

## [1] 0.1742193
[Link]<-[Link]([Link],SAheart$chd,alpha=1)
min([Link]$cvm)

66
## [1] 0.1745402

Pour être plus complet, il faudrait aussi calculer les courbes ROC, ainsi que les tables de prévision sur un
échantillon test. Cette dernière est une information qui peut jouer un rôle. En effet, dans ce problème, on
veut un modèle qui limite au maximum le FPR (“False Positive Rate”), i.e. l’erreur selon laquelle on prévoit
chd = 0 (non atteint par la pathologie) alors que la réalité est chd = 1 (atteint par la pathologie), ce qui est
clairement mieux pour le patient ! Rappelons qu’on peut calculer les prévisions d’un objet de la classe [Link]

avec la commande predict([Link],newx=MatriceTest,s=“[Link]”).

2. On retient le modèle lasso avec interaction. Comme c’est un objet de la classe [Link], on ne peut pas
appliquer la fonction [Link] de la librairie questionr. On procède donc directement (mais sans avoir
de p-valeurs du coup) :

oddLasso<-exp(coef([Link],s="[Link]"))
oddLasso

## 7 x 1 Matrix of class "dgeMatrix"

## s1
## (Intercept) 0.6622834
## tobacco 1.0176694
## ldl 1.0000639
## famhistPresent 0.8491014
## typea 1.0053979
## age 1.0077173
## ldl:famhistPresent 1.0714686

67
SEANCE DE TP6
EXERCICE 1

library(rpart)
library([Link])
head(ptitanic)

## pclass survived sex age sibsp parch

## 1 1st survived female 29.0000 0 0
## 2 1st survived male 0.9167 1 2
## 3 1st died female 2.0000 1 2
## 4 1st died male 30.0000 1 2
## 5 1st died female 25.0000 1 2
## 6 1st survived male 48.0000 0 0

1. Taille du fichier et nombre de survivants.

nrow(ptitanic) #Nb de passagers dans le fichier

## [1] 1309
sum(ptitanic$survived=="survived") #Nb de survivants

## [1] 500

2. Le traitement des “NA” est en dehors du cadre du cours, on les supprime donc. Puis on rebat les données
et on construit les échantillons d’apprentissage et de validation.

ptitanic<-[Link](ptitanic) #Suppression des NA

n<-nrow(ptitanic)
[Link](1234)
perm<-sample(n)
titanic<-ptitanic[perm[1:n],]
taille_app<-(9*n)%/%10
dapp_titanic<-titanic[1:taille_app,] #Echantillon d'apprentissage de taille 90%
dtest_titanic<-titanic[-(1:taille_app),]

3. Arbre maximal :

arbre_max<-rpart(survived~.,data=dapp_titanic,minsplit=5,cp=0)
prp(arbre_max)

68
 yes sex = mal no

pclass = 3rd

age >= 9.5

sibsp >= 2 age >= 57

pclass = 2nd,3rd age >= 1.5

age >= 37 survived survived survived

age >= 32 died sibsp >= 2 pclass = 2nd

age < 32 age >= 18 age >= 25 survived

age < 36 age < 48 age >= 5.5 died

age >= 37 died age >= 55 survived age < 4 age < 3.5

died pclass = 2nd died age < 20 age < 17 parch >= 3

died age < 20 parch >= 1 age < 34 sibsp < 3 survived age < 31 age >= 7

age >= 45 died died died died died died died

age < 45 parch < 1 age < 51 survived survived sibsp < 1 survived

died parch >= 1 died age >= 32 died parch < 1 age < 22

died age >= 15 survived survived parch >= 2 age >= 28 age < 33 parch < 1

sibsp >= 1 died died age >= 50 died died age < 38 age >= 20

age >= 40 survived survived age >= 28 survived parch >= 1 survived age < 15 survived

died died died died died died died died age >= 24

age < 39 age < 26 survived age >= 17 survived

died age < 30 died died died

pclass = 2nd age < 31 survived survived

died died sibsp < 1 sibsp >= 1

survived survived survived survived

died died parch < 1 died

survived

died

Courbe de ses erreurs.

plotcp(arbre_max)

size of tree

1 2 3 5 8 9 12 20 33 41 49 56 64
1.0
X−val Relative Error

0.8
0.6
0.4

Inf 0.14 0.016 0.0048 0.0031 0.0017 5e−04

cp

69
4. Procédure d’élagage de l’arbre maximal.

cp_opti<-arbre_max$cptable[[Link](arbre_max$cptable[,4]),1]
arbre_opti<-prune(arbre_max,cp=cp_opti)
prp(arbre_opti)

yes sex = mal no

age >= 9.5 pclass = 3rd

died sibsp >= 2 died survived

died survived

5. Probabilités de survie/décès des individus de l’échantillon test.

prev<-predict(arbre_opti,newdata=dtest_titanic)
prev

## died survived
## 834 0.82214156 0.17785844
## 211 0.82214156 0.17785844
## 113 0.07142857 0.92857143
## 858 0.82214156 0.17785844
## 479 0.07142857 0.92857143
## 714 0.82214156 0.17785844
## 168 0.07142857 0.92857143
## 884 0.82214156 0.17785844
## 708 0.82214156 0.17785844
## 54 0.82214156 0.17785844
## 912 0.82214156 0.17785844
## 934 0.82214156 0.17785844
## 132 0.07142857 0.92857143
## 1106 0.82214156 0.17785844
## 818 0.82214156 0.17785844
## 690 0.82214156 0.17785844
## 1102 0.93333333 0.06666667
## 429 0.07142857 0.92857143
## 501 0.82214156 0.17785844
## 815 0.82214156 0.17785844
## 1307 0.82214156 0.17785844
## 810 0.82214156 0.17785844
## 50 0.82214156 0.17785844
## 846 0.55714286 0.44285714
## 466 0.07142857 0.92857143

70
## 21 0.82214156 0.17785844
## 550 0.12000000 0.88000000
## 1297 0.82214156 0.17785844
## 349 0.82214156 0.17785844
## 262 0.82214156 0.17785844
## 48 0.82214156 0.17785844
## 722 0.82214156 0.17785844
## 555 0.82214156 0.17785844
## 950 0.82214156 0.17785844
## 552 0.07142857 0.92857143
## 27 0.82214156 0.17785844
## 1245 0.55714286 0.44285714
## 19 0.07142857 0.92857143
## 1270 0.82214156 0.17785844
## 1264 0.82214156 0.17785844
## 98 0.07142857 0.92857143
## 340 0.93333333 0.06666667
## 513 0.82214156 0.17785844
## 777 0.82214156 0.17785844
## 1022 0.82214156 0.17785844
## 77 0.07142857 0.92857143
## 178 0.82214156 0.17785844
## 987 0.82214156 0.17785844
## 842 0.55714286 0.44285714
## 861 0.55714286 0.44285714
## 327 0.82214156 0.17785844
## 567 0.82214156 0.17785844
## 437 0.07142857 0.92857143
## 645 0.82214156 0.17785844
## 913 0.82214156 0.17785844
## 334 0.07142857 0.92857143
## 441 0.07142857 0.92857143
## 427 0.82214156 0.17785844
## 871 0.55714286 0.44285714
## 111 0.82214156 0.17785844
## 412 0.07142857 0.92857143
## 9 0.07142857 0.92857143
## 17 0.82214156 0.17785844
## 1184 0.82214156 0.17785844
## 241 0.82214156 0.17785844
## 521 0.82214156 0.17785844
## 338 0.07142857 0.92857143
## 761 0.82214156 0.17785844
## 393 0.07142857 0.92857143
## 894 0.82214156 0.17785844
## 297 0.07142857 0.92857143
## 1051 0.55714286 0.44285714
## 519 0.82214156 0.17785844
## 568 0.82214156 0.17785844
## 638 0.82214156 0.17785844
## 522 0.07142857 0.92857143
## 630 0.82214156 0.17785844
## 733 0.82214156 0.17785844
## 67 0.07142857 0.92857143

71
## 145 0.07142857 0.92857143
## 783 0.82214156 0.17785844
## 310 0.07142857 0.92857143
## 468 0.07142857 0.92857143
## 196 0.07142857 0.92857143
## 672 0.82214156 0.17785844
## 750 0.55714286 0.44285714
## 658 0.55714286 0.44285714
## 1259 0.55714286 0.44285714
## 746 0.55714286 0.44285714
## 1120 0.82214156 0.17785844
## 788 0.82214156 0.17785844
## 1008 0.55714286 0.44285714
## 791 0.82214156 0.17785844
## 82 0.82214156 0.17785844
## 1229 0.82214156 0.17785844
## 355 0.82214156 0.17785844
## 190 0.82214156 0.17785844
## 854 0.82214156 0.17785844
## 34 0.07142857 0.92857143
## 204 0.82214156 0.17785844
## 798 0.82214156 0.17785844
## 786 0.55714286 0.44285714
## 372 0.07142857 0.92857143
## 237 0.82214156 0.17785844
## 847 0.82214156 0.17785844

On peut comparer ces probabilités à la réalité, par exemple pour le dernier individu de l’échantillon test :

dtest_titanic[n-taille_app,]

## pclass survived sex age sibsp parch

## 847 3rd died male 20 0 0

On peut aussi obtenir la classe prédite en ajoutant le paramètre “type=class”.

class_prev<-predict(arbre_opti,newdata=dtest_titanic,type="class")
class_prev

## 834 211 113 858 479 714 168 884

## died died survived died survived died survived died
## 708 54 912 934 132 1106 818 690
## died died died died survived died died died
## 1102 429 501 815 1307 810 50 846
## died survived died died died died died died
## 466 21 550 1297 349 262 48 722
## survived died survived died died died died died
## 555 950 552 27 1245 19 1270 1264
## died died survived died died survived died died
## 98 340 513 777 1022 77 178 987
## survived died died died died survived died died
## 842 861 327 567 437 645 913 334
## died died died died survived died died survived

72
## 441 427 871 111 412 9 17 1184
## survived died died died survived survived died died
## 241 521 338 761 393 894 297 1051
## died died survived died survived died survived died
## 519 568 638 522 630 733 67 145
## died died died survived died died survived survived
## 783 310 468 196 672 750 658 1259
## died survived survived survived died died died died
## 746 1120 788 1008 791 82 1229 355
## died died died died died died died died
## 190 854 34 204 798 786 372 237
## died died survived died died died survived died
## 847
## died
## Levels: died survived

6. Pour vérifier la pertinence du modèle, il peut aussi être intéressant de dresser la table de confusion sur
l’échantillon test.

table(obs=dtest_titanic$survived,prev=class_prev)

## prev
## obs died survived
## died 55 1
## survived 23 26

Enfin, on peut aussi calculer la probabilité de mauvais classement.

mean(dtest_titanic$survived!=class_prev)

## [1] 0.2285714

7. On compare maintenant ce résultat avec celui d’une régression logistique (donc sans arbre).

library(glmnet)
logit<-glm(survived~.,dapp_titanic,family=binomial)
summary(logit)

##
## Call:
## glm(formula = survived ~ ., family = binomial, data = dapp_titanic)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6741 -0.6957 -0.4227 0.6847 2.5359
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.846766 0.378067 10.175 < 2e-16 ***
## pclass2nd -1.383494 0.242039 -5.716 1.09e-08 ***
## pclass3rd -2.348444 0.240811 -9.752 < 2e-16 ***
## sexmale -2.453698 0.180324 -13.607 < 2e-16 ***

73
## age -0.040431 0.006908 -5.852 4.85e-09 ***
## sibsp -0.374911 0.109311 -3.430 0.000604 ***
## parch 0.078062 0.103327 0.755 0.449961
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1267.89 on 940 degrees of freedom
## Residual deviance: 883.37 on 934 degrees of freedom
## AIC: 897.37
##
## Number of Fisher Scoring iterations: 4

On calcule ensuite le taux de mauvais classement dans l’échantillon test.

prev_logit<-predict(logit,newdata=dtest_titanic,type="response",se=T)
#se=T : pour pouvoir récupérer les probabilités de survie/décès

prevision_logit <- ifelse(prev_logit$fit > 0.5, "survived","died")

#On converti les probabilités de survie/décès en "survived" ou "died"

mean(prevision_logit!=dtest_titanic$survived)

## [1] 0.1904762

74
SEANCE DE TP7

EXERCICE 1
library(kernlab)

##
## Attaching package: 'kernlab'
## The following object is masked from 'package:ggplot2':
##
## alpha
data(spam)
attach(spam)

## The following object is masked _by_ .GlobalEnv:

##
## data
#head(spam)
sum(spam$type=="spam")

## [1] 1813
nrow(spam)

## [1] 4601

1. Commencons avec un modèle de type régression logistique.

logit<-glm(type~.,data=spam,family=binomial)

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

library(boot)
library(glmnet)

[Link](spam,logit,K=10)$delta[1]

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

75
## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred
## [1] 0.05840814

Le warning “[Link]: des probabilités ont été ajustées numériquement à 0 ou 1” n’est pas un bon signe.
Lorsque la log-vraisemblance a été calculée, beaucoup de valeurs égales à 0 ou 1 ont été trouvées. Or, ce
modèle n’est pas fait pour de telles valeurs, mais pour des valeurs strictement comprises entre 0 et 1. Cela
peut engendrer un phénomène d’instabilité dans les estimations. Voyons si on peut quand même améliorer

les prédicteurs précédents en ajoutant, en tout début d’étude, une sélection de variables pas à pas. Faisons-le
uniquement avec un descente pas à pas et le critère BIC, pour aller chercher de la parcimonie (mais le critère
BIC utilisant le log-vraisemblance du modèle, on retrouvera forcément les “warning” qu’on a rencontré plus
haut).

#[Link]<-step(logit,direction="backward",k=log(nrow(spam)))

Partant de ce modèle, on peut ainsi construire le modèle logistique dont on étudie les propriétés :

[Link]<-glm(type~make + address + num3d + our + over + remove + internet +

order + mail + will + addresses + free + business + you +
credit + your + font + num000 + money + hp + hpl + george +
num650 + lab + data + num85 + technology + parts + pm + cs +
meeting + original + project + re + edu + table + conference +
charSemicolon + charExclamation + charDollar + charHash +
capitalLong + capitalTotal,data=spam,family="binomial")

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

library(boot)
library(glmnet)

[Link](spam,[Link],K=10)$delta[1]

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

## Warning: [Link]: fitted probabilities numerically 0 or 1 occurred

76
## [1] 0.05724459

L’amélioration du modèle logistique n’est pas flagrante, et surtout encore quelques warning !

2. Une solution, partant de ce modèle de régression logistique, peut être de pénaliser sa vraisemblance (comme
on l’a vu en pénalisant l’erreur quadratique avec le modèle linéaire) avec une pénalisation de type ridge ou
lasso (pour tenter de rendre le modèle parcimonieux). La syntaxe est la même que pour le modèle linéaire,
à ceci près qu’il faut rajouter les argument family=“binomial” pour faire référence au modèle logistique et
type=“class” pour se ramener à un problème de discrimination. Examinons tout d’abord le modèle ridge.

library(glmnet)
[Link]<-[Link](type~.,data=spam)
ridge<-[Link]([Link],spam$type,family="binomial",alpha=0,[Link]="class")

min(ridge$cvm)

## [1] 0.08758965

Essayons maintenant avec le lasso.

library(glmnet)
lasso<-[Link]([Link],spam$type,family="binomial",alpha=1,[Link]="class")

sum(coef(lasso,s=lasso$[Link])==0) #Nombre de coefficients nuls

## [1] 4
min(lasso$cvm)

## [1] 0.07150619

3. Comparons maintenant avec une forêt d’arbres :

library(randomForest)

## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
##
## margin
foret<-randomForest(type~.,data=spam)
foret

##
## Call:
## randomForest(formula = type ~ ., data = spam)
## Type of random forest: classification
## Number of trees: 500

77
## No. of variables tried at each split: 7
##
## OOB estimate of error rate: 4.5%
## Confusion matrix:
## nonspam spam [Link]
## nonspam 2714 74 0.02654232
## spam 133 1680 0.07335907

La valeur de l’erreur OOB ne permet pas de comparer aux autres prédicteurs. Calculons donc l’erreur K-fold
sur cette forêt :

library(caret)
cv<-trainControl(method="cv",number=10)
rf<-train(type~.,data=spam,method="rf",trControl=cv,tuneGrid=[Link](mtry=round(sqrt(ncol(spam)))))
1-rf$results[rownames(rf$bestTune),]$Accuracy

## [1] 0.04716358
# ou rf$results pour des infos plus complètes

La forêt d’arbre remporte le match !

4. On peut essayer d’améliorer la forêt d’arbres en jouant sur le paramètre mtry (autour de la valeur par
défaut de la fonction randomForest). Compte-tenu du grand nombre de calculs qui seraient nécessaires, on
parallélise le calcul avec la librairie doParallel.

library(caret)
library(doParallel) ## pour paralléliser

## Loading required package: foreach

## Loading required package: iterators
## Loading required package: parallel
[Link] <- [Link](mtry=seq(5,10,by=1))
cv<-trainControl(method="cv",number=10)

registerDoParallel(makePSOCKcluster(4))

[Link] <- train(type~.,data=spam,method="rf",trControl=cv,

tuneGrid=[Link])

[Link](stopCluster(makePSOCKcluster(4)))

[Link]$results

## mtry Accuracy Kappa AccuracySD KappaSD

## 1 5 0.9539163 0.9028612 0.010435783 0.02216272
## 2 6 0.9539163 0.9028854 0.011981743 0.02536667
## 3 7 0.9536984 0.9024611 0.010298289 0.02189884
## 4 8 0.9534815 0.9020229 0.011464656 0.02429363
## 5 9 0.9534829 0.9020780 0.010646580 0.02247289
## 6 10 0.9539172 0.9030663 0.009803885 0.02068778

78
5. La valeur par défaut pour mtry donne donc (presque) le meilleur résultat. On peut aussi reprendre le
même calcul, mais cette fois avec les variables sélectionnées par l’AIC, pour constater qu’il n’y a toujours pas
d’amélioration (et que la valeur mtry par défaut donne -presque- le meilleur prédicteur).

6. On reste donc sur la forêt d’arbres initiale. Etudions maintenant les variables explicatives importantes,
tout d’abord de manière numérique :

foret$importance

## MeanDecreaseGini
## make 8.0263712
## address 12.0717458
## all 25.1390918
## num3d 1.8111581
## our 62.3431631
## over 18.1456957
## remove 175.9677745
## internet 26.1031719
## order 9.2576681
## mail 17.0939929
## receive 20.4876490
## will 23.5707227
## people 8.1152736
## report 4.7713156
## addresses 2.8310559
## free 136.7074220
## business 26.6028802
## email 18.6416924
## you 58.9607872
## credit 11.1494441
## your 129.9928228
## font 5.4083001
## num000 58.3515815
## money 79.3719995
## hp 95.1554531
## hpl 40.5093702
## george 43.9370023
## num650 13.6213695
## lab 5.2414687
## labs 7.6770700
## telnet 4.9131135
## num857 1.4982297
## data 5.9381749
## num415 1.6806936
## num85 7.6338303
## technology 7.6192754
## num1999 29.0797193
## parts 0.9790439
## pm 6.9672805
## direct 2.6182099
## cs 2.3102780
## meeting 11.6273048
## original 3.2268100
## project 4.1021197

79
## re 20.8021405
## edu 40.7041820
## table 0.5028313
## conference 2.7845217
## charSemicolon 11.4436900
## charRoundbracket 25.5534508
## charSquarebracket 6.3821038
## charExclamation 251.9282481
## charDollar 198.5972189
## charHash 7.2528958
## capitalAve 146.9719688
## capitalLong 118.6303199
## capitalTotal 87.5148886

Puis, avec un graphique :

varImpPlot(foret)

foret

charExclamation
charDollar
remove
capitalAve
free
your
capitalLong
hp
capitalTotal
money
our
you
num000
george
edu
hpl
num1999
business
internet
charRoundbracket
all
will
re
receive
email
over
mail
num650
address
meeting

0 50 100 150 200 250

MeanDecreaseGini

80
EXERCICE 2
bank<-[Link]("~/Dropbox/ABIDJAN/DONNEES/[Link]",header=T,sep=" ")
head(bank)

## age job marital education default housing contact

## 14281 36 admin. single [Link] unknown yes cellular
## 32177 39 services married [Link] unknown yes cellular
## 13177 36 admin. married [Link] no yes cellular
## 18265 33 unemployed married [Link] unknown yes cellular
## 36457 55 retired married basic.4y no no cellular
## 29319 42 services married [Link] no no cellular
## day_of_week duration campaign pdays previous poutcome y
## 14281 mon 393 5 999 0 nonexistent no
## 32177 fri 19 6 999 0 nonexistent no
## 13177 wed 570 2 999 0 nonexistent yes
## 18265 wed 590 9 999 0 nonexistent no
## 36457 tue 130 2 999 0 nonexistent no
## 29319 fri 277 2 999 0 nonexistent no
bank$y<-factor(bank$y)

Pas vraiment d’amélioration, mais on retient qu’on peut quand même utiliser le modèle le plus parcimonieux.
Regardons maintenant les modèles ridge et lasso.

library(glmnet)
[Link]<-[Link](y~.,data=bank)
[Link]<-[Link]([Link],bank$y,family="binomial",alpha=0,[Link]="class")

min([Link]$cvm)

## [1] 0.0925
[Link]<-[Link]([Link],bank$y,family="binomial",alpha=1,[Link]="class")

sum(coef([Link],s=lasso$[Link])==0)

## [1] 2
min([Link]$cvm)

## [1] 0.0917

Jusqu’à présent, on conserve donc le modèle logistique. Etudions maintenant un modèle de forêt d’arbres.

library(randomForest)
[Link]<-randomForest(y~.,data=bank)
[Link]

##
## Call:
## randomForest(formula = y ~ ., data = bank)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##

81
## OOB estimate of error rate: 9.27%
## Confusion matrix:
## no yes [Link]
## no 8661 202 0.02279138
## yes 725 412 0.63764292

La valeur de l’erreur OOB ne permet pas de comparer aux autres prédicteurs. On peut néanmoins observer
que le taux de faux positifs est très élevé ! Calculons l’erreur K-fold sur cette forêt :

library(caret)

cv<-trainControl(method="cv",number=10)

[Link]<-train(y~.,data=bank,method="rf",trControl=cv,tuneGrid=[Link](mtry=round(sqrt(ncol(bank)))))

[Link]$results

## mtry Accuracy Kappa AccuracySD KappaSD

## 1 4 0.9064001 0.3512107 0.003860173 0.04215884
Pour les variables explicatives importantes :
varImpPlot([Link])

[Link]

duration
age
job
pdays
day_of_week
education
campaign
poutcome
previous
marital
housing
contact
default

0 100 200 300 400 500 600

MeanDecreaseGini

82