POL6022 - Méthodes quantitatives avancées

Evelyne Brie

Hiver 2026

Apprentissage automatique

Ce lab applique deux méthodes vues en cours sur des données concrètes : la validation croisée pour évaluer un modèle OLS, et les forêts d’arbres décisionnels pour prédire un vote politique.

Packages nécessaires : caret, rpart, rpart.plot, randomForest, carData, class, ggplot2

1. Validation croisée

1.1 Le jeu de données et le modèle

Nous repartons du jeu de données swiss déjà utilisé lors du lab sur les diagnostics OLS. L’objectif change : au lieu d’interpréter les coefficients, on veut évaluer la capacité prédictive du modèle sur de nouvelles données.

data(swiss)
head(swiss)

##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6

Modèle estimé

Nous reprenons le même modèle que lors du lab sur les diagnostics :

\[\widehat{\text{Fertility}}_i = \hat{\beta}_0 + \hat{\beta}_1\,\text{Agriculture}_i + \hat{\beta}_2\,\text{Education}_i + \hat{\beta}_3\,\text{Catholic}_i + \varepsilon_i\]

La question n’est plus “les coefficients sont-ils significatifs?” mais “ce modèle prédit-il bien la fécondité dans des provinces qu’il n’a pas vues?”

1.2 Partition entraînement / test

On divise les données en deux sous-ensembles : 70% pour entraîner le modèle, 30% pour évaluer sa performance sur des observations non vues.

set.seed(123)
n <- nrow(swiss)

# on tire aleatoirement 70% des indices de lignes pour le training set
idx_train <- sample(1:n, size = floor(0.7 * n))

# pour être claire : idx_train n'est qu'une liste de numéros de lignes
idx_train

##  [1] 31 15 14  3 42 37 45 25 26 27  5 38 28  9 29 44  8 39  7 10 34 19  4 41 17
## [26] 11 33 12 46 35 13 21

# les lignes selectionnées vont dans train, les autres dans test
train <- swiss[idx_train, ]
test  <- swiss[-idx_train, ]

nrow(train)

## [1] 32

nrow(test)

## [1] 15

On estime le modèle uniquement sur le training set, puis on prédit sur les deux ensembles.

mod_train <- lm(Fertility ~ Agriculture + Education + Catholic, data = train)

summary(mod_train)

## 
## Call:
## lm(formula = Fertility ~ Agriculture + Education + Catholic, 
##     data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5836  -5.9094   0.1257   3.6321  13.2717 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 87.60468    5.76707  15.191 4.76e-15 ***
## Agriculture -0.20388    0.08599  -2.371   0.0249 *  
## Education   -1.06148    0.17155  -6.188 1.11e-06 ***
## Catholic     0.15790    0.03323   4.752 5.47e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.148 on 28 degrees of freedom
## Multiple R-squared:  0.7292, Adjusted R-squared:  0.7002 
## F-statistic: 25.13 on 3 and 28 DF,  p-value: 4.272e-08

# on prédit sur les deux ensembles avec le même modèle
pred_train <- predict(mod_train, newdata = train)
pred_test  <- predict(mod_train, newdata = test)

# regarder la prédiction pour les 10 premières observations de chaque jeu de données
# (histoire de vous illustrer comment la "prédiction" fonctionne concrètement)
pred_train[1:10]

##      Conthey     Cossonay     Avenches Franches-Mnt    Neuchatel       Sierre 
##     83.71242     68.61354     63.19076     88.95095     52.72056     82.87650 
## V. De Geneve         Oron      Payerne Paysd'enhaut 
##     37.78687     72.40578     68.09311     71.87797

pred_test[1:10]

## Courtelary   Delemont Porrentruy  Echallens   Lausanne     Lavaux     Moudon 
##   72.97357   82.25244   87.27821   74.50107   55.84001   63.61641   73.90006 
##      Nyone       Orbe    Yverdon 
##   66.87990   70.86896   69.98386

# MSE = moyenne des carrés des erreurs de prédiction
mse_train <- mean((train$Fertility - pred_train)^2)
mse_test  <- mean((test$Fertility  - pred_test)^2)

mse_train

## [1] 44.7071

mse_test

## [1] 88.88155

D’autres métriques permettent de mieux comprendre la qualité des prédictions sur le test set.

# Erreur absolue moyenne : dans les mêmes unités que la VD, plus facile à interpréter
mae_test <- mean(abs(test$Fertility - pred_test))
mae_test

## [1] 8.04138

# RMSE : dans les mêmes unités que la VD, mais pénalise davantage les grandes erreurs
rmse_test <- sqrt(mean((test$Fertility - pred_test)^2))
rmse_test

## [1] 9.427701

# Corrélation entre les valeurs prédites et les valeurs réelles dans le test set
# une corrélation proche de 1 indique que le modèle capte bien la structure des données
cor(pred_test, test$Fertility)

## [1] 0.7193623

Interprétation : Le MSE de test est plus élevé que le MSE d’entraînement, ce qui est normal : le modèle a été construit sur le training set, donc il y performe artificiellement mieux. La MAE est plus intuitive que le MSE parce qu’elle est dans les mêmes unités que la variable dépendante : ici, elle indique de combien le modèle se trompe en moyenne sur le taux de fécondité. La corrélation entre valeurs prédites et valeurs réelles mesure à quel point le modèle capture la structure des données : une corrélation élevée indique que les provinces avec une fécondité élevée sont bien classées comme telles par le modèle, même si les valeurs absolues ne sont pas parfaites.

1.3 Validation croisée k-fold

La partition 70/30 dépend du tirage aléatoire : avec un autre set.seed, on obtiendrait des MSE différents. La validation croisée à \(k\) plis résout ce problème en répétant l’évaluation \(k\) fois sur des sous-ensembles différents.

La validation croisée à \(k\) plis répète l’évaluation \(k\) fois sur des sous-ensembles différents. Le choix de \(k\) dépend principalement de la taille de l’échantillon :

\(k = 10\) : valeur standard dans la littérature, bon compromis pour la plupart des tailles d’échantillon
\(k = 5\) : préférable quand \(n\) est petit (moins de 50 observations), comme ici avec swiss
LOOCV (leave-one-out cross-validation) : chaque observation sert tour à tour de test set, à réserver aux très petits échantillons (\(n < 30\))

La règle pratique est que chaque fold de test devrait contenir au moins 5 observations. Avec \(n = 47\), un \(k = 10\) donne environ 5 observations par fold, ce qui est à la limite acceptable. Un \(k = 5\) donnerait environ 9 observations par fold et serait plus stable.

dim(swiss)

## [1] 47  6

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

mod_cv <- train(Fertility ~ Agriculture + Education + Catholic,
                data      = swiss,
                method    = "lm",
                trControl = ctrl)

# RMSE et R2 moyens sur les 10 plis
mod_cv$results

##   intercept     RMSE  Rsquared     MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 7.933204 0.6431886 7.06116 2.574569  0.2804069 2.459848

# Maintenant avec 5 plis
ctrl <- trainControl(method = "cv", number = 20)

mod_cv <- train(Fertility ~ Agriculture + Education + Catholic,
                data      = swiss,
                method    = "lm",
                trControl = ctrl)

# RMSE et R2 moyens sur les 5 plis
mod_cv$results

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 7.369735 0.8278462 6.755252 3.223157  0.3402647 3.083701

Interprétation : Le RMSE de validation croisée est une estimation plus fiable de l’erreur de généralisation que le RMSE d’entraînement, parce qu’il est calculé sur des observations que le modèle n’a jamais vues lors de chaque pli. Le \(R^2\) de CV est un peu plus bas que le \(R^2\) calculé sur les données d’entraînement, ce qui reflète la différence entre fit et généralisation.

Exercice

Ajoutez la variable Examination (pourcentage de recrues ayant obtenu la note la plus élevée à l’examen militaire) au modèle. Comparez le RMSE de validation croisée des deux modèles. Le modèle avec Examination prédit-il mieux?

# Votre code ici

mod_cv2 <- train(Fertility ~ Agriculture + Education + Catholic + Examination,
                 data      = swiss,
                 method    = "lm",
                 trControl = ctrl)

# comparer les RMSE de validation croisée des deux modèles
mod_cv$results$RMSE

## [1] 7.369735

mod_cv2$results$RMSE

## [1] 7.639528

Le RMSE de validation croisée est légèrement plus bas sans Examination, ce qui suggère une meilleure capacité prédictive. En prédiction pure, on privilégie le modèle avec le RMSE de CV le plus bas.

2. Forêts d’arbres décisionnels (facultatif)

2.1 Le jeu de données

Nous utilisons le jeu de données Chile du package carData, qui contient des données d’un sondage politique réalisé en 1988 avant le plébiscite chilien sur la continuation du régime de Pinochet.

data(Chile)
head(Chile)

##   region population sex age education income statusquo vote
## 1      N     175000   M  65         P  35000   1.00820    Y
## 2      N     175000   M  29        PS   7500  -1.29617    N
## 3      N     175000   F  38         P  15000   1.23072    Y
## 4      N     175000   F  49         P  35000  -1.03163    N
## 5      N     175000   F  23         S  35000  -1.10496    N
## 6      N     175000   F  28         P   7500  -1.04685    N

Les variables pertinentes sont :

vote : intention de vote (Y = oui, N = non, A = abstention, U = indécis) – variable dépendante
age : âge du répondant
income : revenu mensuel
statusquo : degré de soutien au statu quo (échelle continue)
sex : sexe (M/F)
education : niveau d’éducation (P = primaire, S = secondaire, PS = post-secondaire)

On se concentre sur les répondants ayant exprimé une intention de vote claire (Y ou N) et on retire les valeurs manquantes.

# on garde seulement les repondants ayant vote Y ou N
chile_bin <- Chile[Chile$vote %in% c("Y", "N"), 
                   c("vote", "age", "income", "statusquo", "sex", "education")]

# droplevels retire les niveaux inutilises (A et U) du facteur vote
chile_bin$vote <- droplevels(chile_bin$vote)
chile_bin      <- na.omit(chile_bin)

nrow(chile_bin)

## [1] 1703

table(chile_bin$vote)

## 
##   N   Y 
## 867 836

On partitionne les données en training (70%) et test (30%).

set.seed(456)
# createDataPartition stratifie le tirage pour conserver la proportion Y/N dans chaque ensemble
idx_chile   <- createDataPartition(chile_bin$vote, p = 0.7, list = FALSE)
chile_train <- chile_bin[idx_chile, ]
chile_test  <- chile_bin[-idx_chile, ]

2.2 Un arbre de décision simple

Avant la forêt, on estime un seul arbre de décision avec rpart pour comprendre la logique de base.

arbre <- rpart(vote ~ age + income + statusquo + sex + education,
               data   = chile_train,
               method = "class")

rpart.plot(arbre, 
           type  = 4, 
           extra = 104,
           main  = "Arbre de decision -- vote chilien")

Interprétation : Chaque noeud indique la condition de division, la classe majoritaire prédite et la proportion d’observations dans ce noeud. La variable statusquo domine l’arbre : le soutien au statu quo est le prédicteur le plus puissant de l’intention de vote. C’est cohérent substantivement : ceux qui appuient le régime en place votent majoritairement Oui.

On évalue les prédictions sur le test set.

pred_arbre <- predict(arbre, newdata = chile_test, type = "class")
confusionMatrix(pred_arbre, chile_test$vote)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   N   Y
##          N 244  30
##          Y  16 220
##                                           
##                Accuracy : 0.9098          
##                  95% CI : (0.8815, 0.9332)
##     No Information Rate : 0.5098          
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.8193          
##                                           
##  Mcnemar's Test P-Value : 0.05527         
##                                           
##             Sensitivity : 0.9385          
##             Specificity : 0.8800          
##          Pos Pred Value : 0.8905          
##          Neg Pred Value : 0.9322          
##              Prevalence : 0.5098          
##          Detection Rate : 0.4784          
##    Detection Prevalence : 0.5373          
##       Balanced Accuracy : 0.9092          
##                                           
##        'Positive' Class : N               
##

2.3 La forêt aléatoire

Un seul arbre est instable : un léger changement dans les données d’entraînement peut donner un arbre très différent. La forêt aléatoire entraîne plusieurs centaines d’arbres sur des sous-échantillons bootstrap et agrège leurs prédictions.

set.seed(456)
foret <- randomForest(vote ~ age + income + statusquo + sex + education,
                      data     = chile_train,
                      ntree    = 500,
                      mtry     = 2,
                      importance = TRUE)
foret

## 
## Call:
##  randomForest(formula = vote ~ age + income + statusquo + sex +      education, data = chile_train, ntree = 500, mtry = 2, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 6.96%
## Confusion matrix:
##     N   Y class.error
## N 566  41  0.06754530
## Y  42 544  0.07167235

Interprétation : L’erreur OOB (Out-Of-Bag) est l’erreur de classification estimée sur les observations non incluses dans chaque bootstrap. C’est une estimation fiable de l’erreur de généralisation sans avoir besoin d’un test set séparé : chaque observation sert de test pour les arbres qui ne l’ont pas vue.

On évalue quand même sur le test set pour rester cohérent avec la section 1.

pred_foret <- predict(foret, newdata = chile_test)
confusionMatrix(pred_foret, chile_test$vote)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   N   Y
##          N 237  26
##          Y  23 224
##                                          
##                Accuracy : 0.9039         
##                  95% CI : (0.875, 0.9281)
##     No Information Rate : 0.5098         
##     P-Value [Acc > NIR] : <2e-16         
##                                          
##                   Kappa : 0.8077         
##                                          
##  Mcnemar's Test P-Value : 0.7751         
##                                          
##             Sensitivity : 0.9115         
##             Specificity : 0.8960         
##          Pos Pred Value : 0.9011         
##          Neg Pred Value : 0.9069         
##              Prevalence : 0.5098         
##          Detection Rate : 0.4647         
##    Detection Prevalence : 0.5157         
##       Balanced Accuracy : 0.9038         
##                                          
##        'Positive' Class : N              
##

Interprétation : Comparez l’exactitude de la forêt avec celle de l’arbre simple. La forêt performe généralement mieux parce qu’elle moyenne les prédictions de centaines d’arbres, ce qui réduit la variance.

2.4 Importance des variables

La forêt aléatoire permet de mesurer la contribution de chaque variable à la performance prédictive.

varImpPlot(foret,
           main = "Importance des variables -- vote chilien")

Interprétation : Le graphique de gauche (Mean Decrease Accuracy) montre de combien l’erreur augmente si on permute aléatoirement les valeurs d’une variable. Le graphique de droite (Mean Decrease Gini) mesure la contribution de chaque variable aux divisions dans les arbres. statusquo domine les deux mesures, ce qui confirme le résultat de l’arbre simple. Rappel : l’importance des variables ne s’interprète pas causalement.

Exercice

Réestimez la forêt avec ntree = 100 et mtry = 4. Comparez l’erreur OOB et l’exactitude sur le test set avec le modèle original (ntree = 500, mtry = 2). Qu’est-ce que cela vous apprend sur le rôle de ces hyperparamètres?

# Votre code ici

set.seed(456)
foret2 <- randomForest(vote ~ age + income + statusquo + sex + education,
                       data       = chile_train,
                       ntree      = 100,
                       mtry       = 4,
                       importance = TRUE)

# erreur OOB a la derniere iteration pour chaque modele
foret$err.rate[500, 1]

##        OOB 
## 0.06957251

foret2$err.rate[100, 1]

##        OOB 
## 0.07544007

pred_foret2 <- predict(foret2, newdata = chile_test)

# exactitude sur le test set pour chaque modèle
mean(pred_foret  == chile_test$vote)

## [1] 0.9039216

mean(pred_foret2 == chile_test$vote)

## [1] 0.8941176

Avec ntree = 100, la forêt est moins stable (les résultats varient plus selon le set.seed). Avec mtry = 4, les arbres sont plus corrélés entre eux car ils choisissent parmi plus de variables à chaque division, ce qui réduit le bénéfice de la diversification. En pratique, on choisit ces hyperparamètres par validation croisée plutôt qu’à la main.