Dans cette section, nous appliquerons la méthode CART pour prédire les délais de paiement en utilisant le langage statistique R et l’IDE RStudio. Nous construirons un arbre de régression sur le jeu de données ModeleDuree pour modéliser la variable continue “DP”. Les données seront divisées en ensembles d’entraînement et de test, l’arbre sera construit avec l’ensemble d’entraînement, et ses performances seront évaluées sur l’ensemble de test.

library(rpart)
library(caret)
library(vip)
library(rpart.plot)
library(partykit)
library(ISLR)
library(dplyr)
library(lubridate)
library(ggparty)
library(readr)
library(ineq)

Introduction

Les méthodes basées sur les arbres segmentent l’espace des prédicteurs en régions simples pour la classification et la régression. Les arbres de décision sont non paramétriques et visent à maximiser l’exactitude en classification ou à minimiser l’erreur en régression. Chaque nœud est défini par une règle “si/alors” basée sur une variable de division et un point de division. Les modèles CART sont préférés aux LM, GLM ou GAM lorsque les relations entre les variables prédictives sont complexes et non linéaires, ou lorsque les interactions entre les variables sont non additives ou nombreuses. Par exemple, nous utiliserons l’ensemble de données Carseats pour illustrer un arbre de classification et un arbre de régression.

Chargement des données.

set.seed(454564)
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.

ModeleDuree <- 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
  )

L’étude de la base de données de prédictions a déjà été réalisée et peut être consultée dans le fichier Analyse des données. Avant de commencer l’analyse, il est nécessaire de séparer la base de données en deux ensembles : une base de test et une base d’apprentissage. À cette fin, la base de test comprendra 30 % des données, tandis que la base d’apprentissage sera constituée des 70 % restants.

Train = sample(1:nrow(ModeleDuree), nrow(ModeleDuree)*0.7, replace = FALSE)

ModeleCARTTrain = ModeleDuree[Train,]
ModeleCARTTest = ModeleDuree[-Train,]

Il est essentiel de régler soigneusement les hyperparamètres avant de lancer le modèle. Une fois cette étape complétée, la modélisation peut commencer..

Réglage des hyperparamètres

Terminologie et éléments de sortie des arbres de décision de cart sont :

  • Nœud racine : Représente l’ensemble de la population d’échantillon et est ensuite divisé en deux ou plusieurs groupes homogènes.

  • split : Le critère utilisé pour diviser un nœud en deux ou plusieurs sous-nœuds.

  • n : Le nombre d’observations dans un nœud.

  • loss : C’est le nombre total de lignes qui seront mal classifiées si la classe prédite pour le nœud est appliquée à toutes les lignes.

  • yval : La prédiction globale pour la branche (Oui ou Non). En général, il s’agit de la valeur de réponse moyenne pour ce sous-ensemble.

  • yprob : La fraction d’observations dans cette branche qui prennent les valeurs Oui et Non.

  • Élagage (Pruning) : Le processus de suppression des sous-nœuds d’un nœud de décision. L’élagage est une solution au problème de surajustement.

  • Complexité : Le paramètre de complexité est utilisé pour établir un niveau de contrôle qui détermine si une division contribue à un meilleur ajustement du modèle. Toute division qui améliore l’ajustement du modèle par un facteur supérieur au paramètre de complexité défini est tentée. La meilleure valeur de cp a la plus faible erreur de validation croisée (xerror).

La complexité cp est automatiquement ajustée dans la fonction rpart à l’aide d’algorithmes. Plusieurs valeurs de cp sont testées, et il suffit de sélectionner celle qui minimise l’erreur xerror.

On commence par construire un arbre avec les paramètres par défaut de rpart

fit0_cart <- rpart(DP~CG+AB+SB+RB+CB+TC+TP+PS+RO, data = ModeleCARTTrain, 
                   parms=list(split="information"), method = "anova")

et on étudie la suite d’arbres emboîtés

printcp(fit0_cart)
## 
## Regression tree:
## rpart(formula = DP ~ CG + AB + SB + RB + CB + TC + TP + PS + 
##     RO, data = ModeleCARTTrain, method = "anova", parms = list(split = "information"))
## 
## Variables actually used in tree construction:
## [1] CG TP
## 
## Root node error: 9695210/291230 = 33.291
## 
## n= 291230 
## 
##         CP nsplit rel error  xerror       xstd
## 1 0.745000      0   1.00000 1.00001 0.00157552
## 2 0.011788      1   0.25500 0.25500 0.00085199
## 3 0.010000      2   0.24321 0.24322 0.00081596

L’erreur d’ajustement (rel.error) diminue normalement, tout comme l’erreur de prévision (xerror). Cela suggère que fit0_cart est suffisamment profond pour être élagué. Par précaution, il peut être utile de considérer un arbre légèrement plus profond.

fit0_cart <- rpart(DP~CG+AB+SB+RB+CB+TC+TP+PS+RO, data = ModeleCARTTrain,
                   minsplit=10,cp=1e-5, method = "anova")
plotcp(fit0_cart)

La valeur de prudence considérée ne donne pas les meilleurs résultats ; au contraire, elle conduit à une erreur relativement élevée.

cptable <- data.frame(fit0_cart$cptable)
subset(cptable, CP == 1e-5)
##        CP nsplit rel.error    xerror         xstd
## 527 1e-05    826 0.1920386 0.2029038 0.0007240722

La valeur choisie donne nsplit = 5346, ce qui entraîne une trop grande division avec un xerror de 0.5387785. La valeur minimale de xerror est atteinte pour cp = 3.707713e-05.

subset(cptable, xerror == min(xerror))
##               CP nsplit rel.error    xerror         xstd
## 305 1.455363e-05    432 0.1967011 0.2016001 0.0007095772

Le modèle avec l’erreur xerror la plus faible est obtenu avec nsplit = 315.

subset(cptable, rel.error == min(rel.error))
##        CP nsplit rel.error    xerror         xstd
## 527 1e-05    826 0.1920386 0.2029038 0.0007240722

Cependant, pour l’erreur rel.error, notre valeur de cp de prudence minimise l’erreur avec un très grand nombre de divisions. Nous conservons donc la valeur de cp qui offre la plus petite xerror.

Regression sur les données

cp.opt <- subset(cptable, xerror == min(xerror))$CP

On recupere le modele ajusté avec la meilleure valeur de cp.

fit0_cart.opt <- prune(fit0_cart, cp=cp.opt)

Visualisation des arbres de regression

rpart.plot(fit0_cart.opt)

Il est difficile d’examiner l’arbre de régression en raison de son grand nombre de branches. Faire un zoom sur les 10 premières branches permet d’observer :

prp(fit0_cart.opt, type = 3, cex = 0.8, extra = 101, main = "Arbre de Régression (Limité)")

Il est également difficile d’observer les variables de découpage de l’arbre.

#subtree <- subset(fit_party, depth <= 3) 

#ggparty(subtree) +
#  theme_minimal() +
#  coord_cartesian(clip = 'off')

Importance des variables

La visualisation de l’arbre aide à interpréter l’algorithme mais ne hiérarchise pas les variables par importance. Une mesure, le score d’importance, évalue chaque variable selon son utilité dans l’arbre. Chaque coupure maximise le gain d’impureté. Le score d’importance d’une variable est la somme des gains d’impureté qu’elle génère.

Soit \(T\) un arbre avec \(|T|\) nœuds terminaux et \(|T| - 1\) nœuds internes. Le gain d’impureté au nœud interne \(t\) est noté \(g_t\), et \(X_c(t)\) est la variable utilisée pour cette coupure. L’importance de la variable \(X_j\) est :

\[ I_j(T) = \sum_{t=1}^{|T|-1} g_t \cdot 1_{c(t)=j} \]

La fonction rpart calcule automatiquement ces scores pour l’arbre optimal.

vip(fit0_cart.opt)

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 grand groupe de code (CG).

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(fit0_cart.opt, ModeleCARTTest)
actuals <- ModeleCARTTest$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.8033088

Le \(R^2\) de 0.8033088 indique que notre modèle de GLM n’explique qu’environ 80.33% de la variance totale des données. Cela suggère que le modèle capture une bonne partie de la relation entre les variables explicatives et la variable dépendante. C’est é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.588005

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

Mean Relative Error (MRE)

le Mean Relative Error (MRE) 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, par valeur réelles. Il est défini par :

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

Le MRE 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.

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

MRE(actuals, predictions)
## [1] 0.8391236

Un MRE 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] 1.988904

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.0009656672

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 <- 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(fit0_cart.opt, newdata = ModeleCARTTrain) - ModeleCARTTrain$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] 3.544497e-15
var(residuals)
## [1] 6.548314

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.