# Pour les tracés
library(ggplot2)

1 Utilisation de données de test

1.1 Chargement du jeu de données

On reprend ci-après le jeu de données déjà utilisé plusieurs fois dans les exemples.

load(url("https://irma.math.unistra.fr/~jberard/donnees_swautoins_ext_sim"))
str(donnees)
## 'data.frame':    239444 obs. of  7 variables:
##  $ Kilometres: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ Zone      : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ Bonus     : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Make      : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Claims    : int  0 1 3 3 6 1 1 1 2 2 ...
##  $ Insured   : num  6.36 7.37 6.81 6.56 13.53 ...
##  $ Payment   : num  0 5475 21866 15230 31354 ...
donnees[["Zone"]] <- relevel(donnees[["Zone"]], ref = "4")
donnees[["Bonus"]] <- relevel(donnees[["Bonus"]], ref = "7")
donnees[["Make"]] <- relevel(donnees[["Make"]], ref = "7-9")

1.2 Séparation en données d’ajustement et de test

On décide maintenant de diviser le jeu de données en deux parties : un jeu de données d’ajustement, comportant \(80\%\) des données, et un jeu de données de test, comportant \(20\%\) des données. La division est effectuée en utilisant un tirage (pseudo-)aléatoire uniforme parmi les sous-ensembles possibles.

# Nombre de lignes du jeu de données d'ajustement
nb_ajust <- ceiling(0.8 * nrow(donnees)) 
# Fixation de la graine du générateur pseudo-aléatoire pour assurer la reproductibilité
set.seed(123456) 
# Tirage aléatoire des indices pour l'ajustement
indices_ajust <- sample(x = (1:nrow(donnees)), size = nb_ajust) 
# Extraction des données utilisées pour l'ajustement
donnees_ajust <- donnees[indices_ajust,] 
# Extraction des données utilisées pour le test
donnees_test <- donnees[- indices_ajust,] 
# Vérification de la proportion d'exposition ainsi obtenue
sum(donnees_ajust[["Insured"]]) / sum(donnees[["Insured"]])
## [1] 0.7999452

1.3 Ajustement du modèle sur les données d’ajustement

On procède ensuite à l’ajustement d’un modèle GLM log-Poisson pour la fréquence de sinistres en utilisant un terme d’offset pour prendre en compte les différences d’exposition. Le modèle est ajusté en faisant uniquement appel aux données d’ajustement uniquement.

ajustement <- glm(Claims ~ Kilometres + Zone + Bonus + Make + offset(log(Insured)),
                  family = poisson(link = "log"),
                  data = donnees_ajust)

1.4 Mesures de performance prédictive sur les données de test

On commence par la déviance.

# Extraction de la fonction de déviance
fonc_deviance <- ajustement$family$dev.resids
# Calcul des espérances modélisées sur les données de test
mu_pred_test <- predict(ajustement, newdata = donnees_test, type = "response")
# Extraction des réponses observées sur les données de test
reponse_test <- donnees_test[["Claims"]]
# Normalisation par l'exposition des valeurs précédentes
expo_test <- donnees_test[["Insured"]]
mu_pred_test_n <- mu_pred_test / expo_test
reponse_test_n <- reponse_test / expo_test
# Calcul des déviances sur chaque ligne du jeu de données de test
dev_test <- fonc_deviance(reponse_test_n, mu_pred_test_n, wt = expo_test)
# Déviance totale sur les données de test
dev_test_tot <- sum(dev_test)
dev_test_tot
## [1] 44347.7
# Déviance moyenne sur les données de test
dev_test_moy <- mean(dev_test) 
dev_test_moy
## [1] 0.9260712

On continue avec l’erreur quadratique.

# Erreur quadratique sur chaque ligne du jeu de données de test
eq_test <- expo_test * (reponse_test_n -  mu_pred_test_n)^2
# Erreur quadratique totale sur les données de test
eq_test_tot <- sum(eq_test)
eq_test_tot
## [1] 2291.419
# Erreur quadratique moyenne sur les données de test
eq_test_moy <- mean(eq_test)
eq_test_moy
## [1] 0.04784954

1.5 Comparaison avec les valeurs moyennes obtenues sur les données d’ajustement

On calcule maintenant les mêmes quantités qu’auparavant, mais sur les données d’ajustement. (On reprend ligne-à-ligne le code précédent pour souligner le parallélisme entre les deux, mais, dans le cas des données d’ajustement, on peut bien entendu obtenir ces valeurs de manière plus directe.)

On commence avec la déviance.

mu_pred_ajust <- predict(ajustement, newdata = donnees_ajust, type = "response")
reponse_ajust <- donnees_ajust[["Claims"]]
expo_ajust <- donnees_ajust[["Insured"]]
mu_pred_ajust_n <- mu_pred_ajust / expo_ajust
reponse_ajust_n <- reponse_ajust / expo_ajust
dev_ajust <- fonc_deviance(reponse_ajust_n, mu_pred_ajust_n, wt = expo_ajust)
dev_ajust_tot <- sum(dev_ajust)
dev_ajust_tot
## [1] 176289.6
dev_ajust_moy <- mean(dev_ajust) 
dev_ajust_moy
## [1] 0.920303

Puis l’erreur quadratique.

eq_ajust <- expo_ajust * (reponse_ajust_n -  mu_pred_ajust_n)^2
# Erreur quadratique totale sur les données de test
eq_ajust_tot <- sum(eq_ajust)
eq_ajust_tot
## [1] 9120.283
# Erreur quadratique moyenne sur les données de test
eq_ajust_moy <- mean(eq_ajust)
eq_ajust_moy
## [1] 0.04761158

On constate que, pour les quantités calculées, les valeurs moyennes obtenues sur les données de test et sur les données d’ajustement sont très voisines, ce qui suggère que l’on n’est pas ici dans une situation de sur-ajustement du modèle. On note que les valeurs obtenues sur les données d’ajustement sont à chaque fois (très légèrement) inférieures aux valeurs obtenues sur les données de test, ce qui est cohérent.

De manière précise, les écarts relatifs entre les valeurs obtenues sont les suivants :

(dev_test_moy - dev_ajust_moy) / dev_ajust_moy
## [1] 0.006267741
(eq_test_moy - eq_ajust_moy) / eq_ajust_moy
## [1] 0.004998014

1.6 Comparaison avec les résultats obtenus pour un autre modèle

A titre d’illustration, on calcule les mêmes mesures de performance prédictive pour un modèle GLM de même forme mais n’utilisant pas la variable Bonus.

# Comparaison avec un modèle n'utilisant pas la variable Bonus
ajustement_sans_Bonus <- glm(Claims ~ Kilometres + Zone + Make + offset(log(Insured)), 
                             family = poisson(link = "log"),
                             data = donnees_ajust)
mu_pred_test_sans_Bonus <- predict(ajustement_sans_Bonus, newdata = donnees_test, type = "response")
mu_pred_test_n_sans_Bonus <- mu_pred_test_sans_Bonus / expo_test
dev_test_sans_Bonus <- fonc_deviance(reponse_test_n, mu_pred_test_n_sans_Bonus, wt = expo_test)
dev_test_tot_sans_Bonus <- sum(dev_test_sans_Bonus)
dev_test_tot_sans_Bonus
## [1] 48754.01
# Déviance moyenne sur les données de test
dev_test_moy_sans_Bonus <- mean(dev_test_sans_Bonus) 
dev_test_moy_sans_Bonus
## [1] 1.018084

Pour rappel, les valeurs obtenues pour le modèle incluant la variable Bonus étaient les suivantes :

dev_test_tot
## [1] 44347.7
dev_test_moy
## [1] 0.9260712

On constate ainsi la nette dégradation des performances prédictives lorsque la variable Bonus est omise.

On reprend le calcul avec l’erreur quadratique.

eq_test_sans_Bonus <- expo_test*(reponse_test_n -  mu_pred_test_n_sans_Bonus)^2
# Erreur quadratique totale sur les données de test
eq_test_tot_sans_Bonus <- sum(eq_test_sans_Bonus)
eq_test_tot_sans_Bonus
## [1] 2563.854
# Erreur quadratique moyenne sur les données de test
eq_test_moy_sans_Bonus <- mean(eq_test_sans_Bonus)
eq_test_moy_sans_Bonus
## [1] 0.05353854

Les valeurs obtenues pour le modèle incluant la variable Bonus étaient les suivantes :

eq_test_tot
## [1] 2291.419
eq_test_moy
## [1] 0.04784954

2 Autres critères : validation croisée et AIC

ajustement <- glm(Claims ~ Kilometres + Zone + Bonus + Make + offset(log(Insured)),
                  family = poisson(link = "log"),
                  data = donnees)

Calculer de manière explicite le critère de validation croisée “leave-one-out” représenterait ici un temps de calcul prohibitif, car cela supposerait de réajuster le modèle autant de fois qu’il y a de lignes (\(239444\)).

# Validation croisée leave-one-out pour la déviance
resultat_loocv <- outilscoursglm::cv.glm_modif(data = donnees, glmfit = ajustement, K = nrow(donnees))

On préfèrera donc utiliser, par exemple, l’approximation basée sur les leviers.

# Validation croisée pour la déviance en utilisant l'approximation par les leviers
resultat_approx_cv <- outilscoursglm::approx_loocv_deviance_leviers(ajustement)
resultat_approx_cv
## [1] 220679

On calcule maintenant un critère de validation croisée “K-fold” avec \(K=10\).

# Validation croisée K-fold avec K=10
resultat_10_fold_cv <- outilscoursglm::cv.glm_modif(data = donnees, glmfit = ajustement, K = 10)
resultat_10_fold_cv
## [1] 220685.6

On calcule maintenant le critère de validation croisée généralisée.

# Validation croisée généralisée GCV*
resultat_gcv <- outilscoursglm::glm_deviance_GCV(ajustement)
resultat_gcv
## [1] 220679.1

Et pour finir le critère AIC, d’abord extrait par une fonction dédiée, puis recalculé manuellement.

# AIC
resultat_AIC <- AIC(ajustement)
resultat_AIC
## [1] 404488.2
# Calcul manuel
- 2 * as.numeric(logLik(ajustement)) + 2 * ajustement$rank
## [1] 404488.2

A titre d’illustration, on calcule de nouveau les quantités précédentes, pour un modèle n’incluant pas la variable Bonus.

ajustement_sans_Bonus <- glm(Claims ~ Kilometres + Zone + Make + offset(log(Insured)),
                             family = poisson(link = "log"),
                             data = donnees)
resultat_approx_cv_sans_Bonus <- outilscoursglm::approx_loocv_deviance_leviers(ajustement_sans_Bonus)
resultat_approx_cv_sans_Bonus
## [1] 243094.3
resultat_10_fold_cv_sans_Bonus <- outilscoursglm::cv.glm_modif(data = donnees, glmfit = ajustement_sans_Bonus, K = 10)
resultat_10_fold_cv_sans_Bonus
## [1] 243089.9
resultat_gcv_sans_Bonus <- outilscoursglm::glm_deviance_GCV(ajustement_sans_Bonus)
resultat_gcv_sans_Bonus
## [1] 243094.4
resultat_AIC_sans_Bonus <- AIC(ajustement_sans_Bonus)
resultat_AIC_sans_Bonus
## [1] 426899.4

Pour rappel, les valeurs obtenues pour le modèle incluant la variable Bonus étaient les suivantes :

resultat_approx_cv
## [1] 220679
resultat_10_fold_cv
## [1] 220685.6
resultat_gcv
## [1] 220679.1
resultat_AIC
## [1] 404488.2

La comparaison est sans appel, et montre, pour chacune des mesures utilisées, la nette dégradation des performances prédictives lorsque la variable Bonus est omise.

3 Courbes de Lorenz

On reprend la séparation en données d’ajustement et données de test utilisée au début de ce bloc-notes, en procédant de nouveau à l’ajustement du modèle sur les données d’ajustement.

ajustement <- glm(Claims ~ Kilometres + Zone + Bonus + Make + offset(log(Insured)), 
                family = poisson(link = "log"),
                data = donnees_ajust)

On procède maintenant au calcul et au tracé de la courbe de Lorenz associée, calculée sur les données de test.

# Calcul de la courbe de Lorenz 
courbe_Lorenz <- outilscoursglm::concentration_curve_glm(ajustement = ajustement,
                                                         donnees = donnees_test,
                                                         nom_reponse = "Claims",
                                                         nom_exposition = "Insured",
                                                         type_rep_modele = "brut")
# Tracé
ggplot() + 
  geom_line(data = courbe_Lorenz, aes(x = cumul_exposition, y = cumul_reponse)) +
  theme_bw() +
  labs(title = "Courbe de Lorenz sur donn\u{00E9}es test", 
       x = "Exposition cumul\u{00E9}e",
       y = "R\u{00E9}ponse cumul\u{00E9}e") +
  geom_line(data = data.frame(x = c(0,1), y = c(0,1)), aes(x = x, y = y), col = "red", linetype = 2) +
  xlim(0,1) +
  ylim(0,1)

A titre d’illustration, on trace également la courbe obtenue avec un modèle n’utilisant pas la variable Bonus.

# Comparaison avec un modèle n'utilisant pas la variable Bonus
ajustement_sans_Bonus <- glm(Claims ~ Kilometres + Zone + Make + offset(log(Insured)),
                             family = poisson(link = "log"),
                             data = donnees_ajust)
courbe_Lorenz_2 <- outilscoursglm::concentration_curve_glm(ajustement = ajustement_sans_Bonus,
                                                           donnees = donnees_test,
                                                           nom_reponse = "Claims",
                                                           nom_exposition = "Insured",
                                                           type_rep_modele = "brut")
ggplot() + 
  geom_line(data = courbe_Lorenz, aes(x = cumul_exposition, y = cumul_reponse, colour = "complet")) +
  geom_line(data = courbe_Lorenz_2, aes(x = cumul_exposition, y = cumul_reponse, colour = "sans Bonus")) +
  theme_bw() +
  labs(title = "Courbe de Lorenz sur donn\u{00E9}es test", 
       x = "Exposition cumul\u{00E9}e",
       y = "R\u{00E9}ponse cumul\u{00E9}e") +
  geom_line(data = data.frame(x = c(0,1), y = c(0,1)), aes(x = x, y = y), col = "red", linetype = 2) +
  xlim(0,1) +
  ylim(0,1) +
  scale_colour_manual(name = "", values = c("complet" = "black", "sans Bonus" = "blue"))

On ajoute également un tracé de courbe de Lorenz “extrême”, à titre de comparaison.

# Courbe de Lorenz extrême (réponse moyenne prédite = réponse observée)
courbe_Lorenz_extr <- outilscoursglm::concentration_curve_devin(donnees = donnees_test,
                                                               nom_reponse = "Claims",
                                                               nom_exposition = "Insured")
ggplot() + 
  geom_line(data = courbe_Lorenz, aes(x = cumul_exposition, y = cumul_reponse, colour = "complet")) +
  geom_line(data = courbe_Lorenz_2, aes(x = cumul_exposition, y = cumul_reponse, colour = "sans Bonus")) +
  geom_line(data = courbe_Lorenz_extr, aes(x = cumul_exposition, y = cumul_reponse, colour = "clairvoyant")) +
  theme_bw() +
  labs(title = "Courbe de Lorenz sur donn\u{00E9}es test", 
       x = "Exposition cumul\u{00E9}e",
       y = "R\u{00E9}ponse cumul\u{00E9}e") +
  geom_line(data = data.frame(x = c(0,1), y = c(0,1)), aes(x = x, y = y), col = "red", linetype = 2) +
  xlim(0,1) +
  ylim(0,1) +
  scale_colour_manual(name = "", values = c("complet" = "black", "sans Bonus" = "blue", "clairvoyant" = "green"))

On effectue maintenant le calcul des aires délimitées par les différentes courbes, ainsi que les rapports entre celles-ci.

# Calculs d'aires
aire_sous_courbe <- outilscoursglm::aire_sous_courbe(courbe_Lorenz)
aire_sous_courbe
## [1] 0.3637649
aire_sous_courbe_2 <- outilscoursglm::aire_sous_courbe(courbe_Lorenz_2)
aire_sous_courbe_2
## [1] 0.4303665
aire_sous_courbe_extr <- outilscoursglm::aire_sous_courbe(courbe_Lorenz_extr)
aire_sous_courbe_extr
## [1] 0.1361847
ratio_aire <- (0.5 - aire_sous_courbe) / (0.5 - aire_sous_courbe_extr)
ratio_aire
## [1] 0.3744622
ratio_aire_2 <- (0.5 - aire_sous_courbe_2) / (0.5 - aire_sous_courbe_extr)
ratio_aire_2
## [1] 0.1913978

4 Tracé de courbe ROC en régression logistique

On illustre ci-après le tracé de courbe ROC en régression logistique, sur un exemple tiré de [David Hosmer, Stanley Lemeshow, and Rodney Sturdivant, Applied logistic regression, John Wiley & Sons, 2013.]

On commence par charger le jeu de données d’ajustement et le jeu de données de test.

load(url("https://irma.math.unistra.fr/~jberard/BURN1000")) 
load(url("https://irma.math.unistra.fr/~jberard/BURN_EVAL_1")) 
str(donnees)
## 'data.frame':    1000 obs. of  11 variables:
##  $ ID      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ FACILITY: int  11 1 12 1 1 6 22 1 1 1 ...
##  $ DEATH   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ AGE     : num  26.6 2 22 37.3 52.1 50.2 2.5 53.8 31.9 41.1 ...
##  $ GENDER  : int  1 0 0 1 1 1 0 0 1 1 ...
##  $ RACEC   : int  1 0 0 1 1 1 0 1 1 1 ...
##  $ TBSA    : num  25.3 5 2 2 6 7 7 0.9 2 22 ...
##  $ INH_INJ : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ FLAME   : int  1 0 0 0 1 0 0 1 0 1 ...
##  $ TBSAFP1 : num  5.03 2.24 1.41 1.41 2.45 ...
##  $ AGEFP1  : num  7.08 0.04 4.84 13.91 27.14 ...
str(donnees_test)
## 'data.frame':    500 obs. of  11 variables:
##  $ ID      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ FACILITY: int  15 48 62 32 28 28 55 15 20 15 ...
##  $ DEATH   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ AGE     : num  26 48.8 15.8 38.2 0.5 37.1 24.3 17.1 37.5 15 ...
##  $ GENDER  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ RACEC   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ TBSA    : num  10 3 4 8 1 1 8 10 0.9 10.3 ...
##  $ INH_INJ : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ FLAME   : int  0 1 0 1 0 0 1 1 0 1 ...
##  $ TBSAFP1 : num  3.16 1.73 2 2.83 1 ...
##  $ AGEFP1  : num  6.76 23.8144 2.4964 14.5924 0.0025 ...

On ajuste ensuite le modèle de régression logistique proposé, puis l’on calcule les espérances (ici, ce sont des probabilités) prédites par le modèle sur les données de test.

# Régression logistique binaire
ajustement <- glm(DEATH ~ 
                    AGEFP1 + TBSAFP1 + RACEC + INH_INJ,
                  family = binomial(link = "logit"),
                  data = donnees) 
# Valeurs moyennes prédites sur les données de test
val_pred_test <- predict(ajustement, newdata = donnees_test, type = "response")

On procède maintenant au tracé de la courbe ROC (en calculant également l’aire sous la courbe).

# Paquet pour la manipulation des courbes ROC
library(pROC) 
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
# Calcul de la courbe ROC sur les données de test
courbe_ROC <- roc(response = donnees_test[["DEATH"]], predictor = val_pred_test) 
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Tracé de la courbe ROC
plot(courbe_ROC, grid = TRUE, lwd = 1, identity.lty = 2, identity.col = "red") 

# Aire sous la courbe ROC
auc(courbe_ROC) 
## Area under the curve: 0.9663