Dans ce document, nous mettrons en place un modèle GLM (Generalized Linear Model) sur nos données pour modéliser et prédire les délais de paiement. Nous allons également implémenter des modèles de régression régularisée, tels que Lasso, Ridge et Elastic Net, afin d’améliorer la performance du modèle. Pour cette étude, nous utiliserons les bibliothèques suivantes :

library(readr)
library(ineq)  # Pour la fonction de la courbe de Lorenz
library(haven)
library(MASS)
library(glmnet)
library(ggplot2)
library(Metrics)
library(mgcv)
library(FactoMineR)
library(caret)
library(dplyr)

Presentation des bases et codage des variables.

La base à utiliser est la suivante :

MODELEDURE <- read_csv("C:/Users/JKLEVINE/Desktop/Memoire/MODELEDURE.csv")
names(MODELEDURE) <- c("CG","AB","SB","RB","CB","TC","TP","PS","PP","DP","RO")
head(MODELEDURE, 5)
## # A tibble: 5 × 11
##   CG        AB SB    RB    CB    TC    TP       PS    PP    DP    RO
##   <chr>  <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 ANHOSP     0 F     BRE   E     IN    I        10    11    32  29.9
## 2 ANHOSP     0 F     BRE   E     IN    V        10    11    31  26.5
## 3 ANHOSP     0 F     CVL   E     CO    I         6     7    48  67.5
## 4 ANHOSP     0 F     CVL   E     CO    V         6     9   105  26.5
## 5 ANHOSP     0 F     HDF   E     CO    I        10    11    45  48.2

Les modalités sont codées en nombres entiers en fonction de leur fréquence totale dans la base de données, comme mentionné dans le document Analyse des données. Par exemple, pour la variable du sexe, nous aurions :

  • Sexe F : 1
  • Sexe M : 0

De même pour les catégories de bénéficiaires :

  • Catégorie Conjoint : 0
  • Catégorie Enfant : 1
  • Catégorie Assuré : 2

Etc.

On obtient la nouvelle base de données suivante :

MODELEDURE <- MODELEDURE %>%
  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
  )
ModeleGLM <- MODELEDURE # copie de la base de données

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.

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

ModeleGLMTrain = ModeleGLM[Train,]
ModeleGLMTest = ModeleGLM[-Train,]

Une fois cette separation faite, la modelisation peut commencer.

Regression sur les données

fit0_GLM = glm(DP~CG+AB+SB+RB+CB+TC+TP+PS+RO, data = ModeleGLMTrain)

Lors de l’analyse statistique, la variable age du bénéficiaire AB présente des p-valeurs supérieures à α = 5%, ce qui indique qu’il n’y a pas de preuve statistiquement significative suggérant que ces variables ont un effet différent de zéro sur la réponse étudiée.

summary(fit0_GLM)
## 
## Call:
## glm(formula = DP ~ CG + AB + SB + RB + CB + TC + TP + PS + RO, 
##     data = ModeleGLMTrain)
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  5.513e+00  1.043e-02  528.798   <2e-16 ***
## CG           4.656e-01  1.723e-03  270.178   <2e-16 ***
## AB          -1.745e-06  1.196e-04   -0.015    0.988    
## SB          -5.132e-01  4.671e-03 -109.856   <2e-16 ***
## RB          -1.839e-01  6.321e-04 -291.021   <2e-16 ***
## CB          -2.402e-01  4.139e-03  -58.040   <2e-16 ***
## TC          -2.752e-01  4.921e-03  -55.930   <2e-16 ***
## TP           1.051e+01  4.812e-03 2183.820   <2e-16 ***
## PS          -1.433e-02  2.214e-04  -64.711   <2e-16 ***
## RO          -1.474e-03  1.476e-05  -99.856   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 7.525223)
## 
##     Null deviance: 48566353  on 1456152  degrees of freedom
## Residual deviance: 10957801  on 1456143  degrees of freedom
## AIC: 7071291
## 
## Number of Fisher Scoring iterations: 2

Un deuxième modèle, sans cette variable, est lancé pour essayer de déterminer si son exclusion améliore la performance du modèle ou si elle a un impact significatif sur les résultats.

fit0_GLM_sans_AB = glm(DP~CG+CB+SB+RB+TC+TP+PS+RO, data = ModeleGLMTrain)
summary(fit0_GLM_sans_AB)
## 
## Call:
## glm(formula = DP ~ CG + CB + SB + RB + TC + TP + PS + RO, data = ModeleGLMTrain)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.513e+00  6.957e-03  792.50   <2e-16 ***
## CG           4.656e-01  1.723e-03  270.18   <2e-16 ***
## CB          -2.402e-01  3.571e-03  -67.25   <2e-16 ***
## SB          -5.132e-01  4.671e-03 -109.87   <2e-16 ***
## RB          -1.839e-01  6.275e-04 -293.15   <2e-16 ***
## TC          -2.752e-01  4.847e-03  -56.79   <2e-16 ***
## TP           1.051e+01  4.804e-03 2187.24   <2e-16 ***
## PS          -1.433e-02  2.214e-04  -64.72   <2e-16 ***
## RO          -1.474e-03  1.475e-05  -99.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 7.525218)
## 
##     Null deviance: 48566353  on 1456152  degrees of freedom
## Residual deviance: 10957801  on 1456144  degrees of freedom
## AIC: 7071289
## 
## Number of Fisher Scoring iterations: 2

La différence au niveau de l’AIC est de 2 ; le deuxième modèle a un AIC inférieur de 1 par rapport au précédent. Cependant, pour déterminer quelles variables doivent être prises en compte dans le modèle, une sélection par AIC (StepAIC) est mise en place avecle tout premier modele.

fit_stepAIC = stepAIC(fit0_GLM, direction = c("both"))
## Start:  AIC=7071291
## DP ~ CG + AB + SB + RB + CB + TC + TP + PS + RO
## 
##        Df Deviance     AIC
## - AB    1 10957801 7071289
## <none>    10957801 7071291
## - TC    1 10981342 7074414
## - CB    1 10983151 7074654
## - PS    1 10989313 7075471
## - RO    1 11032836 7081227
## - SB    1 11048619 7083308
## - CG    1 11507115 7142515
## - RB    1 11595137 7153612
## - TP    1 46846110 9186812
## 
## Step:  AIC=7071289
## DP ~ CG + SB + RB + CB + TC + TP + PS + RO
## 
##        Df Deviance     AIC
## <none>    10957801 7071289
## + AB    1 10957801 7071291
## - TC    1 10982070 7074509
## - PS    1 10989319 7075470
## - CB    1 10991834 7075803
## - RO    1 11032873 7081229
## - SB    1 11048635 7083308
## - CG    1 11507117 7142514
## - RB    1 11604507 7154786
## - TP    1 46958670 9190305

Le meilleur modèle est en effet celui qui ne possède pas la variable AB.

fit1 <- fit_stepAIC

Toutes les variables sauf AB sont importantes car l’exclusion de chacune augmente l’AIC, ce qui suggère que chaque variable contribue de manière significative à la qualité du modèle. Les augmentations de l’AIC varient, ce qui peut indiquer l’importance relative des variables.

Performance du modele

Pour juger de la performance de notre modèle prédictif, il existe plusieurs indicateurs. En posant \(k\) le nombre de paramètres à estimer dans le modèle, \(n\) la taille de l’échantillon, \(L\) la vraisemblance du modèle estimée, \(\mathcal{L}\) sa log-vraisemblance, et \(\mathcal{L}_{\text{sat}}\) la log-vraisemblance du modèle saturé, on a les indicateurs suivants :

  • La déviance : \(D = -2 (\mathcal{L} - \mathcal{L}_{\text{sat}})\)
  • AI : \(\text{AIC} = 2k - 2 \ln(L)\)
  • BIC : \(\text{BIC} = -2\mathcal{L} + \ln(n)k\)

Le critère AIC (Akaike Information Criterion), proposé par Hirotugu Akaike en 1973 , le critère d’information bayésien (BIC) introduit par Gideon Schwarz en 1978 et la déviance sont des indicateurs les plus utilisés dans les modèles GLM. Ils permettent de comparer plusieurs modèles, des valeurs faibles de ces derniers indiquent un bon ajustement.

Une analyse graphique des résidus et des résidus de Pearson est recommandée pour valider le modèle.

fit1$deviance
## [1] 10957801
AIC(fit1)
## [1] 7071289
BIC(fit1)
## [1] 7071411

Le modèle est certes meilleur, mais ses indicateurs restent tout de même assez élevés, ce qui peut être expliqué par la taille des données. D’autres indicateurs permettent de juger des performances prédictives du modèle. les indicateurs qui suivent sont calculé sur la base de données test et validations.

predictions <- predict(fit1, ModeleGLMTest)
actuals <- ModeleGLMTest$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.7766997

Le \(R^2\) de 0.7766997 indique que notre modèle de GLM n’explique qu’environ 77.67% de la variance totale des données. Bien que cela puisse sembler satisfaisant pour un modèle de GLM, cela reste relativement faible 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.745291

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] 1.018142

Un MRE de 58 % reste peu convaincant, car cela signifie que jusqu’à 48 % 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.133153

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

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.

predicted <- predict(fit1)
residuals <- ModeleGLMTrain$DP - predicted
plot(Lc(residuals), main="Courbe de Lorenz des Résidus")

plot(Lc(predicted), 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.

Validaiton croisée

Une autre possibilité est de faire appel à la validation croisée de type « K-fold », qui ne nécessite qu’un nombre limité de ré-estimations du modèle (et est applicable avec tout modèle de régression, indépendamment d’une structure linéaire sous-jacente).

On crée une partition (en général par tirage aléatoire) de \(\{1, \ldots, N\}\) en \(K\) sous-ensembles \(I_1, \ldots, I_K\) de tailles approximativement égales, et l’on considère

\[ \rho_{\text{K-fold}} = \sum_{k=1}^{K} \sum_{i \in I_k} \rho(y_i, \hat{\mu}^{[-k]}_K, w_i), \]

où \(\hat{\mu}^{[-k]}_K\) est la valeur moyenne prédite pour \(Y_i\), par le modèle ajusté uniquement sur les données \((y', x')_{` \in \{1, \ldots, N\} \setminus I_k}\).

On utilise typiquement \(K \in \{5, \ldots, 10\}\). Noter que de nombreuses extensions et…

control <- trainControl(method = "cv", number = 10)  # 10-fold cross-validation
model <- caret::train(DP ~ CG + RB + CB + TC + TP + PS + RO + SB, data = ModeleGLMTrain, 
               method = "glm", family = "gaussian", trControl = control)
print(model$results)
##   parameter     RMSE  Rsquared      MAE      RMSESD  RsquaredSD       MAESD
## 1      none 2.743419 0.7743415 2.132174 0.006955959 0.001132637 0.004151532
summary(model)
## 
## Call:
## NULL
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.513e+00  6.957e-03  792.50   <2e-16 ***
## CG           4.656e-01  1.723e-03  270.18   <2e-16 ***
## RB          -1.839e-01  6.275e-04 -293.15   <2e-16 ***
## CB          -2.402e-01  3.571e-03  -67.25   <2e-16 ***
## TC          -2.752e-01  4.847e-03  -56.79   <2e-16 ***
## TP           1.051e+01  4.804e-03 2187.24   <2e-16 ***
## PS          -1.433e-02  2.214e-04  -64.72   <2e-16 ***
## RO          -1.474e-03  1.475e-05  -99.88   <2e-16 ***
## SB          -5.132e-01  4.671e-03 -109.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 7.525218)
## 
##     Null deviance: 48566353  on 1456152  degrees of freedom
## Residual deviance: 10957801  on 1456144  degrees of freedom
## AIC: 7071289
## 
## Number of Fisher Scoring iterations: 2
print(model$resample)
##        RMSE  Rsquared      MAE Resample
## 1  2.746220 0.7742657 2.130177   Fold01
## 2  2.751606 0.7727512 2.137032   Fold02
## 3  2.733265 0.7757599 2.125034   Fold03
## 4  2.739337 0.7747903 2.129874   Fold04
## 5  2.741891 0.7747889 2.133164   Fold05
## 6  2.733127 0.7761475 2.127481   Fold06
## 7  2.751119 0.7729229 2.136729   Fold07
## 8  2.750451 0.7734070 2.137035   Fold08
## 9  2.740140 0.7747623 2.131285   Fold09
## 10 2.747036 0.7738189 2.133934   Fold10

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 <- fit1$residuals
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] -1.840598e-11
var(residuals)
## [1] 7.525177

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.

Regulation

En régression linéaire, des variables fortement corrélées ou un nombre de variables supérieur aux observations peuvent entraîner des estimations imprécises et un sur-ajustement. L’estimateur des coefficients par les MCO est : \(\hat{\beta} = (X'X)^{-1}X'Y\) mais variables corrélées rendent \(X'X\) non inversible, augmentant le risque d’overfitting. La régularisation contrôle les valeurs des coefficients, les rendant plus stables et réduisant le sur-ajustement.

Regression Ridge

La régression Ridge impose une contrainte sur la norme \(l_2\) des paramètres, en la forçant à être inférieure ou égale à une valeur \(\lambda\). L’objectif est de minimiser :

\[ \sum_{i=1}^{n} [y_i - (\beta_{0} + \beta_{1}x_{i,1} + \beta_{2}x_{i,2} + \dots + \beta_{p}x_{i,p})]^2 \quad \text{ sous la contrainte } \|\beta\|^{2}_{2} \leq \lambda \]

Ce problème équivaut à minimiser la fonction de coût suivante :

\[ \sum_{i=1}^{n} [y_i - (\beta_{0} + \beta_{1}x_{i,1} + \beta_{2}x_{i,2} + \dots + \beta_{p}x_{i,p})]^2 + \lambda \|\beta\|^{2}_{2} \]

où \(\|\beta\|^{2}_{2}\) est la norme \(l_2\) des coefficients \(\beta\), et \(\lambda\) est le paramètre de régularisation.

La solution analytique pour les coefficients Ridge est :

\[ \hat{\beta}_{ridge} = (X'X + \lambda I_p)^{-1}X'Y \]

Bien que l’estimateur Ridge soit biaisé, il présente une variance plus faible comparée à celle de la régression linéaire classique. Il est recommandé de centrer et de réduire les données avant l’application de la régression Ridge pour équilibrer les pénalisations entre les variables. De plus, la régression Ridge tend à donner des coefficients similaires pour les variables fortement corrélées.

X <- ModeleGLMTrain %>%
  dplyr::select(CG, CB, RB, TC, TP, PS, RO, SB) %>%
  as.matrix()

Y <- ModeleGLMTrain %>%
  dplyr::select(DP) %>%
  as.matrix()

X_test <- ModeleGLMTest %>%
  dplyr::select(CG, CB, RB, TC, TP, PS, RO, SB) %>%
  as.matrix()

Y_test <- ModeleGLMTest %>%
  dplyr::select(DP) %>%
  as.matrix()
model_ridge <- glmnet::cv.glmnet(X, Y, alpha = 0)
plot(model_ridge)

(best_lambdaridge <- model_ridge$lambda.min)
## [1] 0.4985629

La valeur du coefficient est relativement élevée, proche de 50 %, ce qui est typique d’un modèle Elastic Net. Le modèle s’éloigne donc du modèle GLM classique.

coefficients <- coef(model_ridge, s = "lambda.min")
print(coefficients)
## 9 x 1 sparse Matrix of class "dgCMatrix"
##                      s1
## (Intercept)  5.86623460
## CG           0.38019653
## CB          -0.20656112
## RB          -0.17449118
## TC          -0.26650708
## TP           9.63872953
## PS          -0.01328816
## RO          -0.00126932
## SB          -0.47144242
predictions <- predict(model_ridge, newx = X, s = "lambda.min")
predictions <- as.vector(predictions)


actual <- Y  # Assurez-vous que Y est également un vecteur si nécessaire
ss_total <- sum((actual - mean(actual))^2)
ss_residual <- sum((actual - predictions)^2)
r_squared <- 1 - (ss_residual / ss_total)


print(paste("R²:", r_squared))
## [1] "R²: 0.769110106891583"

On retrouve une valeur du \(R^2\) assez proche de celle du GLM simple, mais plus faible.

predictions <- as.vector(predict(model_ridge, X_test))
actuals <- as.vector(Y_test)

print(R2(actuals, predictions))
## [1] 0.6556058
print(relative_error(actuals, predictions))
## [1] 0.000715984
print(MRE(actuals, predictions))
## [1] 1.128472
print(MAE(actuals, predictions))
## [1] 2.176373
print(RMSE(actuals, predictions))
## [1] 2.775761

Le \(R^2\) évalué sur les données de test est plus faible que celui évalué sur les données d’entraînement, et donc inférieur à celui du modèle GLM simple. Le modèle Ridge ne performe pas mieux que le modèle GLM. On observe également un sous-apprentissage, comme le montre l’erreur relative positive.

Regression Lasso

La régression lasso est une régression linéaire régularisée utilisant la norme \(l_1\) au lieu de la norme euclidienne (\(l_2\)). Lasso signifie Least Absolute Shrinkage and Selection Operator.

On minimise :

\[ \sum_{i=1}^{n} [y_i - (\beta_{0} + \beta_{1}x_{i,1} + \beta_{2}x_{i,2} + \dots + \beta_{p}x_{i,p})]^2 \quad \text{sous } \|\beta\|_{1} = |\beta_{0}| + |\beta_{1}| + \dots + |\beta_{p}| \leq \lambda.\]

Contrairement à la régression ridge, la régression lasso n’a pas d’expression explicite pour \(\hat{\beta}\) et nécessite un algorithme itératif de descente. Le lasso produit des modèles parcimonieux en mettant certains coefficients à zéro, ce qui permet la sélection de variables et la réduction de dimension.

model_lasso <- glmnet::cv.glmnet(X, Y, alpha = 1)
plot(model_lasso)

(best_lambdaLasso <- model_lasso$lambda.min)
## [1] 0.01074121

Le coefficient \(\lambda\) pour le Lasso est assez faible, ce qui signifie que le modèle Lasso est assez proche d’un modèle GLM classique.

coefficients <- coef(model_lasso, s = "lambda.min")
print(coefficients)
## 9 x 1 sparse Matrix of class "dgCMatrix"
##                       s1
## (Intercept)  5.476113286
## CG           0.455774553
## CB          -0.217436506
## RB          -0.180886874
## TC          -0.246863339
## TP          10.482101335
## PS          -0.013261824
## RO          -0.001388303
## SB          -0.489417046
predictions <- predict(model_lasso, newx = X_test, s = "lambda.min")
predictions <- as.vector(predictions)
actual <- as.vector(Y_test)

On retrouve une valeur du \(R^2\) assez proche de celle du GLM simple.

print(R2(actuals, predictions))
## [1] 0.7720368
print(relative_error(actuals, predictions))
## [1] 0.0006677653
print(MRE(actuals, predictions))
## [1] 1.026823
print(MAE(actuals, predictions))
## [1] 2.134779
print(RMSE(actuals, predictions))
## [1] 2.745431

Le R² basé sur les données de test est meilleur pour ce modèle que pour le modèle Ridge, mais on reste avec un sous-apprentissage car l’erreur relative est positive.

Regression Elastic Net

La régression Elastic Net combine les avantages du Lasso et du Ridge. Le Lasso crée des modèles parcimonieux et interprétables, tandis que le Ridge évite le sur-apprentissage et regroupe les variables corrélées. Elastic Net minimise :

\[ \sum_{i=1}^{n} [y_i - (\beta_{0} + \beta_{1}x_{i,1} + \beta_{2}x_{i,2} + \dots + \beta_{p}x_{i,p})]^2 + \lambda (\alpha \|\beta\|_{1} + \frac{1}{2}(1-\alpha) \|\beta\|^{2}_{2}) \]

avec \(\alpha\) entre [0,1]. Si \(\alpha = 1\), on retrouve le Lasso, et si \(\alpha = 0\), on obtient le Ridge.

model_ElasticNet <- glmnet::cv.glmnet(X, Y, alpha = 0.5)
plot(model_ElasticNet)

(best_lambdaElasticNet <- model_ElasticNet$lambda.min)
## [1] 0.01957399

Ce \(\lambda\) est plus faible.

predictions <- predict(model_ElasticNet, newx = X_test, s = "lambda.min")
predictions <- as.vector(predictions)

Malgré les regulations, les R² est toujours le même.

print(R2(actuals, predictions))
## [1] 0.7697644
print(relative_error(actuals, predictions))
## [1] 0.0006686629
print(MRE(actuals, predictions))
## [1] 1.028363
print(MAE(actuals, predictions))
## [1] 2.135161
print(RMSE(actuals, predictions))
## [1] 2.745438

Ce modele se rapproche plus du modele Lasso.

GAM (Modèle Additif Général)

Nous utilisons un modèle Additif Général (GAM) pour prédire les délais de paiement en utilisant plusieurs variables explicatives en essayant de capté les lien non lineaire.

model_gam <- gam(DP ~ CG + RB + CB + TC + TP + PS + RO + SB, data = ModeleGLMTrain)
predictions <- predict(model_gam, newdata = ModeleGLMTest)

Prediction des valeurs

print(R2(actuals, predictions))
## [1] 0.7766997
print(relative_error(actuals, predictions))
## [1] 0.0006659586
print(MRE(actuals, predictions))
## [1] 1.018142
print(MAE(actuals, predictions))
## [1] 2.133153
print(RMSE(actuals, predictions))
## [1] 2.745291

ce modele reste assez proche de celui du modele GLM classique.