Les forêts aléatoires, introduites par Breiman (2001), agrègent des arbres construits sur des échantillons bootstrap, où les variables de coupure sont sélectionnées parmi un sous-ensemble aléatoire, augmentant ainsi la diversité et réduisant la corrélation entre les arbres. Cette méthode améliore la précision des prédictions et réduit le surapprentissage. Les techniques de machine learning nécessitent le calibrage des paramètres pour optimiser les modèles, et les sections suivantes se concentreront sur cette optimisation pour les forêts aléatoires.

library(readr)
library(caret)
library(ranger)
library(randomForestExplainer)
library(dplyr)
library(kernlab)
library(ineq)
library(doParallel) 
library(randomForest)

Preparation des données

set.seed(1234)
ModeleDuree <- read_csv("C:/Users/JKLEVINE/Desktop/Memoire/MODELEDURE.csv")
names(ModeleDuree) <- c("CG","AB","SB","RB","CB","TC","TP","PS","PP","DP","RO")

Les modalités catégorielles sont codées en valeurs entières selon leur fréquence d’apparition.

ModeleFA <- ModeleDuree %>%
  mutate(
    SB = match(SB, names(sort(table(SB), decreasing = TRUE))) - 1, 
    CB = match(CB, names(sort(table(CB), decreasing = TRUE))) - 1,
    CG = match(CG, names(sort(table(CG), decreasing = TRUE))) - 1,
    TC = match(TC, names(sort(table(TC), decreasing = TRUE))) - 1,
    TP = match(TP, names(sort(table(TP), decreasing = TRUE))) - 1,
    RB = match(RB, names(sort(table(RB), decreasing = TRUE))) - 1
  )

Nous allons tenter de predire des des nouvelles valeurs des delais de paiement en se basant sur celle de cette base de données.

Separation des données en données de test et d’apprentissage

On sépare l’échantillon en un échantillon d’apprentissage de taille 3 000 et un échantillon de validation de taille 1 601. L’échantillon d’apprentissage sera utilisé pour construire la forêt et sélectionner ses paramètres, celui de validation pour estimer les performances de la forêt et vérifier que cette dernière ne fait pas de sur-ajustement.

set.seed(1234)
Train = sample(1:nrow(ModeleFA), nrow(ModeleFA)*0.7, replace = FALSE)

ModeleFATrain = ModeleFA[Train,]
ModeleFATest = ModeleFA[-Train,]

Forêts aléatoires

Tout comme les arbres, les forêts aléatoires permettent de prédire une variable quantitative ou qualitative à partir de variables explicatives quantitatives et qualitatives. Cette famille d’algorithmes a été introduite par Breiman (2001) et permet de pallier, dans une certaine mesure, le manque de stabilité des arbres. La méthode est simple à mettre en œuvre et possède souvent de bonnes performances en termes de qualité de prédiction sur des jeux de données complexes, notamment en présence d’un grand nombre de variables explicatives.

Comme son nom l’indique, une forêt aléatoire consiste à agréger un grand nombre d’arbres de classification ou de régression (selon la nature de la variable à expliquer). Le caractère aléatoire du processus vient tout d’abord du fait que les arbres sont construits sur des échantillons bootstrap. Les échantillons bootstrap s’obtiennent généralement par tirage de \(n\) observations avec remise parmi \(n\) dans l’échantillon initial et les arbres sont construits selon une légère variante de l’algorithme CART présenté. En particulier, un second niveau d’aléatoire est introduit à chaque étape de la construction des arbres : à chaque nœud, seul un sous-ensemble de variables est sélectionné aléatoirement pour choisir la coupure. L’estimateur des forêts aléatoires en régression s’écrit :

\[ F(x) = \frac{1}{B} \sum_{k=1}^{B} T_k(x) \] où \(T_k(x)\) représente l’arbre de décision construit sur le \(k\)-ème échantillon bootstrap.

Algorithme d’une Forêt Aléatoire

On a l’agorithme suivant :

Entrées :

  • B : un entier positif.

  • mtry : un entier entre 1 et d.

  • min.node.size : un entier plus petit que n.

Pour b allant de 1 à B :

  1. Tirer aléatoirement avec remise un échantillon de taille n dans {1, ..., n}. Noter Ib l’ensemble des indices sélectionnés et Dn,b = {(xi, yi), i ∈ Ib} l’échantillon bootstrap associé.

  2. Construire un arbre CART à partir de Dn,b :

    1. Choisir mtry variables au hasard parmi les d variables explicatives.

    2. Sélectionner la meilleure coupure Xj ≤ s parmi les mtry variables sélectionnées.

    3. Ne pas découper un nœud s’il contient moins de min.node.size observations.

  3. Noter T(., θb, Dn) l’arbre obtenu.

Retourner :

\[ f_n(x) = \frac{1}{B} \sum_{b=1}^B T(x, θ_b, D_n) \]

Calibration d’un algorithme avec caret

En apprentissage supervisé, l’objectif est de prédire une sortie \(y \in Y\) à partir d’entrées \(x = (x_1, \ldots, x_d) \in \mathbb{R}^d\) en trouvant un algorithme de prévision \(f : \mathbb{R}^d \to Y\). Les forêts aléatoires sont une méthode puissante pour cette tâche. Elles utilisent un ensemble d’arbres de décision pour améliorer la précision et réduire le sur-ajustement.

Pour chaque algorithme, y compris les forêts aléatoires, il est crucial de calibrer les paramètres, tels que le nombre d’arbres ou la profondeur maximale des arbres. Une calibration appropriée permet d’équilibrer la complexité du modèle et d’éviter le sur-ajustement.

Le package caret facilite cette calibration en offrant une approche générale pour optimiser les paramètres des forêts aléatoires, parmi plus de 230 algorithmes. Nous examinerons comment utiliser la fonction train du package pour sélectionner automatiquement les meilleurs paramètres pour les forêts aléatoires.

Construire et analyser une forêt aléatoire

Plusieurs packages permettent de construire des forêts aléatoires, nous utiliserons le package randomForest et la fonction randomForest. La syntaxe de cette fonction est similaire aux autres fonctions de régression ou de discrimination, une formule indique la variable à expliquer et les variables explicatives.

foretRF <- randomForest(DP~CG+AB+SB+RB+CB+TC+TP+PS+RO,data=ModeleFATrain)
foretRA <- ranger(DP~CG+AB+SB+RB+CB+TC+TP+PS+RO,data=ModeleFATrain)
foretRF
## 
## Call:
##  randomForest(formula = DP ~ CG + AB + SB + RB + CB + TC + TP +      PS + RO, data = ModeleFATrain) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 6.965863
##                     % Var explained: 50.64

La fonction ranger du package ranger est plus rapide que randomForest et offre des performances améliorées pour les grandes ensembles de données et les grands nombres d’arbres. ranger est optimisé pour une exécution rapide et peut traiter des jeux de données plus volumineux avec une efficacité accrue.

foretRA
## Ranger result
## 
## Call:
##  ranger(DP ~ CG + AB + SB + RB + CB + TC + TP + PS + RO, data = ModeleFATrain) 
## 
## Type:                             Regression 
## Number of trees:                  500 
## Sample size:                      14371 
## Number of independent variables:  9 
## Mtry:                             3 
## Target node size:                 5 
## Variable importance mode:         none 
## Splitrule:                        variance 
## OOB prediction error (MSE):       6.965963 
## R squared (OOB):                  0.5063836

La sortie informe que la forêt aléatoire ajustée est une regression avec 500 arbres agrégés et un paramètre mtry de 3. Le taux d’erreur et la matrice de confusion sont estimés par la méthode Out Of Bag (OOB), qui utilise les observations non incluses dans les échantillons bootstrap pour évaluer la performance.

l’erreur de prediction est 8.056835%.

foretRA$prediction.error
## [1] 6.965963

Lors de la construction d’une forêt aléatoire avec randomForest, environ 1/3 des observations sont laissées hors sac (OOB) pour chaque arbre, et oob.times enregistre combien de fois chaque observation a été utilisée comme OOB. Ces observations sont utilisées pour évaluer l’erreur de prédiction de la forêt.

Sélectionner les paramètres de la forêt

Pour construire une forêt aléatoire, il est crucial de calibrer correctement le nombre d’arbres (ntree) et le nombre de variables (mtry) sélectionnées à chaque étape.

Selection du ntree

Un grand nombre d’arbres est recommandé, et la stabilité des taux d’erreur peut être vérifiée en utilisant une fonction de tracé.

plot(foretRF)

La courbe noire indique l’erreur de régression totale. Les taux d’erreur sont stables pour un nombre d’arbres entre 400 et 500, donc une forêt de 500 arbres peut être conservée. Si ce nombre d’arbres est insuffisant, le paramètre ntree peut être ajusté dans la fonction randomForest.

Selection du mtry

Pour le paramètre mtry, nous allons utiliser le package caret en utilisant une grille de candidat.

grille.mtry <- data.frame(mtry=seq(1,15,by=2))

Le taux de bien classés (Accuracy) est estimé en utilisant la technique Out Of Bag (OOB), qui est particulièrement adaptée aux forêts aléatoires et plus rapide qu’une validation croisée. Avec cette méthode, nous choisissons mtry = 10.

ctrl <- trainControl(method="oob")
cl <- makePSOCKcluster(4)
registerDoParallel(cl)
sel.mtry <- caret::train(DP~CG+AB+SB+RB+CB+TC+TP+PS+RO,data=ModeleFATrain,
                         method="rf",trControl=ctrl,
tuneGrid=grille.mtry)
sel.mtry
## Random Forest 
## 
## 14371 samples
##     9 predictor
## 
## No pre-processing
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared 
##    1    2.888453  0.4087526
##    3    2.639604  0.5062396
##    5    2.662181  0.4977572
##    7    2.676108  0.4924885
##    9    2.686621  0.4884931
##   11    2.689510  0.4873927
##   13    2.689387  0.4874394
##   15    2.686315  0.4886098
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 3.

La meilleur valeurs reste celle que nous avions de base, mtry = 3.

Contruction du modele

Une fois les hyperparamètres optimaux identifiés pour améliorer le modèle, la prochaine étape consiste à mettre en œuvre ce modèle et à évaluer sa performance.

ModeleFA <- foretRF

Importance des variables

L’importance des variables dans une forêt aléatoire est calculée en comparant l’erreur OOB avec l’erreur obtenue lorsque les valeurs d’une variable sont permutées aléatoirement. Un écart important indique une grande importance de la variable. Pour visualiser les 10 variables les plus importantes, utilisez la fonction importance et un diagramme en bâtons.

var.imp <- ModeleFA$importance
var.imp.normalized <- 100 * var.imp / sum(var.imp)
ord <- order(var.imp.normalized, decreasing = TRUE)
barplot(var.imp.normalized[ord][1:10],
        names.arg = rownames(var.imp)[ord][1:10],
        cex.names = 0.6,
        ylab = "Importance Normalisée (%)",
        main = "Top 10 Variables les Plus Importantes")

La variable TP se distingue nettement avec le score d’importance le plus élevé. Elle est essentielle pour prédire les délais de paiement. Elle est suivie de loin par le montant de remboursement (RO) et le l’age du beneficiaire (AB).

Le type de paiement (TP) et le type de code acte (CG) jouent un rôle crucial dans la détermination des montants des paiements. Le montant de remboursement (RO) a également un impact significatif sur les délais de paiement, soulignant son importance dans notre modèle de prédiction.

En revanche, le sexe du bénéficiaire SB et le type de contrat TC ont une influence marginale sur les prédictions, montrant qu’ils contribuent peu aux résultats du modèle.

Performance du modele

Nous allons utiliser plusieur indicateurs pour calculer les performances predictives du modele.

predictions <- predict(ModeleFA, ModeleFATest)
actuals <- ModeleFATest$DP

Le \(R^2\)

Le coefficient de détermination \(R^2\) représente le pourcentage de variance expliquée par le modèle. Il est défini par :

\[ R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \]

  • \(R^2 = 1\) signifie que le modèle prédit parfaitement les données : les résidus sont nuls et tous les points sont parfaitement alignés sur la droite (ou la courbe) dont la pente est non nulle.
  • \(R^2 = 0\) signifie que le modèle n’a aucun pouvoir prédictif.

Dans de nombreux cas, \(R^2\) est égal au carré du coefficient de corrélation \(R\).

R2=function(actuals, predictions){
  R2 <- var(predictions)/var(actuals)
  return(R2)
}

print(R2(actuals, predictions))
## [1] 0.5170881

Le \(R^2\) de 0.5307152 indique que notre modèle de GLM n’explique qu’environ 53.07 % de la variance totale des données. Cela suggère que le modèle ne capture qu’a peine plus de la moitiée la relation entre les variables explicatives et la variable dépendante. C’est legerement élévé, mais pas assez suffisant compte tenu de la précision à laquelle nous aspirons.

Erreur quadratique moyenne (RMSE)

Le RMSE (Root Mean Squared Error) est une mesure de l’erreur moyenne entre les valeurs prédites par le modèle et les valeurs réelles de la variable cible. Il calcule la racine carrée de la moyenne des carrés des différences entre les prédictions et les valeurs réelles. Il est défini par :

\[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \]

Le RMSE est souvent utilisé pour évaluer à quel point les prédictions du modèle sont proches des valeurs réelles. Un RMSE plus bas indique une meilleure performance du modèle, car cela signifie que les prédictions sont plus proches des valeurs réelles.

RMSE <- function(actuals, predictions){
  RMSE <- sqrt((mean((actuals - predictions)^2)))
  return(RMSE)
}
                 
print(RMSE(actuals, predictions))         
## [1] 2.606178

Notre RMSE est assez faible, c’est assez bon.

Le MAPE (Mean Absolute Percentage Error)

Le MAPE (Mean Absolute Percentage Error), ou Erreur Absolue Moyenne en Pourcentage, est une mesure de la précision d’un modèle de prédiction. Il calcule la moyenne des erreurs absolues relatives entre les valeurs prédites et les valeurs réelles, exprimées en pourcentage. Il est défini par :

\[ \text{MAPE} = \frac{1}{n} \sum_{i=1}^{n} \frac{\left| y_i - \hat{y}_i \right|}{y_i} \]

Le MAPE est utile pour évaluer la performance d’un modèle en termes de pourcentage d’erreur, ce qui facilite la comparaison entre différents modèles ou différentes séries de données.

MAPE=function(actuals, predictions){
  relative_errors = abs((actuals - predictions)) / actuals
  MAPE = mean(relative_errors)
  return(MAPE)
}

MAPE(actuals, predictions)
## [1] 0.4540684

Un MAPE de 44% % reste peu convaincant, car cela signifie que jusqu’à 44 % des valeurs prédites peuvent contenir des erreurs significatives par rapport aux valeurs réelles. En d’autres termes, près de la moitié des données pourrait être mal prédite, ce qui indique une performance relativement faible du modèle.

MAE (Erreur Absolue Moyenne)

Le MAE, ou Erreur Absolue Moyenne, est une mesure de la précision d’un modèle de prédiction. Il calcule la moyenne des valeurs absolues des erreurs entre les valeurs prédites et les valeurs réelles. L’erreur absolue pour chaque observation est définie comme suit :

\[ \text{MAE} = \frac{ \sum_{i=1}^{N} |y_{i} - \hat{y}_{i}|}{N} \]

Le MAE est la moyenne de ces erreurs absolues sur toutes les observations et donne une indication générale de l’erreur moyenne en unités de la variable mesurée.

MAE=function(actuals, predictions){
  absolute_errors <- abs(actuals - predictions)
  MAE <- mean(absolute_errors)
  return(MAE)
  }

print(MAE(actuals, predictions))
## [1] 2.014981

L’erreur relative

L’erreur relative mesure la différence entre les valeurs prédites et les valeurs réelles en termes proportionnels. Elle est calculée en utilisant la formule suivante :

\[ \text{Erreur relative} = \frac{\sum_{i=1}^{n} \left| y_i - \hat{y}_i \right|}{\sum_{i=1}^{n} y_i} \]

Cette formule calcule la somme des erreurs absolues pour chaque observation, divisée par la somme des valeurs réelles. L’erreur relative donne une mesure proportionnelle des erreurs par rapport aux valeurs réelles.

relative_error <- function(actuals, predictions){
  relative_error <- sum(actuals - predictions)/sum(actuals)
  return(relative_error)
}

print(relative_error(actuals, predictions))
## [1] -0.00839534

L’erreur relative positive est interprétée comme une sous-estimation. Dans ce cas, nous observons donc une sous-estimation.

Courbe de Lorenz

La courbe de Lorenz fournit une visualisation de la capacité du modèle à discriminer efficacement entre les valeurs moyennes plus ou moins élevées de la réponse (dans le cas où les réponses sont dans l’intervalle [0, +∞[). Classiquement, on considère un échantillon de test \((y_i, x_i)_{i \in I_{\text{Test}}}\) indépendant de l’échantillon utilisé pour l’ajustement du modèle. La courbe est obtenue de la manière suivante :

  1. Tri des données de test par ordre croissant de la valeur moyenne (normalisée) prédite pour obtenir l’échantillon trié \((y(i), x(i))_{i \in I_{\text{Test}}}\).
  2. Tracé des points de coordonnées \((u(i), v(i))_{i \in I_{\text{Test}}}\) définies par :
    • \(u(i)\) : somme des expositions associées aux données \((y(\cdot), x(\cdot))\) telles que \(\cdot \leq i\)
    • \(v(i)\) : somme des réponses \((y(\cdot))\) telles que \(\cdot \leq i\)

En général, les coordonnées u et v sont normalisées respectivement par la somme des expositions et des réponses totales.

residuals <- as.vector(actuals - predictions)
plot(Lc(residuals), main="Courbe de Lorenz des Résidus")

plot(Lc(predictions), main="Courbe de Lorenz des Valeurs Prédites")

La courbe de Lorenz (en noir plus foncé) représente la part de patrimoine détenue par les bénéficiaires lorsqu’ils sont classés par ordre croissant de délais de paiement. Plus les courbes s’éloignent de la diagonale (en noir), plus la distribution est inégale. Dans notre cas, les courbes montrent que les distributions ne sont pas particulièrement inégalitaires.

Analyse des Résidus

L’examen graphique des résidus vise à :

  • Vérifier l’adéquation entre le modèle et les données : Cela permet de s’assurer que le modèle est bien adapté aux données observées.
  • Identifier d’éventuelles tendances ou structures globales non prises en compte : L’analyse des résidus peut révéler des motifs ou des structures dans les données qui n’ont pas été capturés par le modèle.
  • Repérer les données « aberrantes » : L’analyse permet de détecter les observations qui s’écartent considérablement des prévisions du modèle, ce qui peut indiquer des anomalies ou des erreurs dans les données.
residuals <- predict(foretRF, newdata = ModeleFATest) - ModeleFATest$DP
hist(residuals, 
     main = "Histogramme des résidus de fit1 avec densité normale", 
     xlab = "Résidus",
     col = "skyblue",
     border = "white",
     prob = TRUE  # Utilisation de la densité de probabilité
)
curve(dnorm(x, mean = mean(residuals), sd = sd(residuals)), 
      add = TRUE, col = "darkblue", lwd = 2)

# 2. QQ-plot
par(mfrow = c(1, 2))  # Pour afficher deux graphiques côte à côte
qqnorm(residuals)
qqline(residuals, col = "red")
#png("C:/Users/JKLEVINE/Desktop/JUDE KLEVINE/qq_plot_residus.png", width = 800, height = 600)
# Réinitialiser le paramètre de par
par(mfrow = c(1, 1))

mean(residuals)
## [1] 0.05982225
var(residuals)
## [1] 6.789686

Les résidus sont alignés sur la droite rouge et semblent suivre une distribution normale. Nous avons une moyenne nulle, mais la variance de notre modèle est différente de 1.