RGP
Attéméné Kouassi
2022-09-04
- INTRODUCTION
- PARTIE I : TRAITEMENT DES DONNEES BRUTES
- PARTIE II : MODELISATION DE LA PROBABILITE DE SURVENANCE DES ACCIDENTS
- CONCLUSION
Package utilisés
library("missForest")
library("DataExplorer")
library("questionr")
library("VIM")
library("dplyr")
library("DescTools")
library("grDevices")
library("zoo")
library("nnet")
library("mfx")
library("rhandsontable")
library("visdat")
library("naniar")
library("rpart")
library("gtsummary")
library("GGally")INTRODUCTION
La régression linéaire multiple est une méthode de régression mathématique pour décrire les variations d’une variable endogène associée aux variations de plusieurs variables exogènes. La désignation “multiple” fait référence au fait qu’il y a plusieurs variables explicatives xj pour expliquer y. La désignation “linéaire” correspond au .. L’analyse par régression linéaire multiple est une des solutions qui existe pour observer les liens entre une variable quantitative dépendante et n variables quantitatives. L’objectif étant d’établir la relation entre Y et les autres variables.
Notre étude portera sur la modélisation du cout des sinistres accidents pour une compagnie d’assurances. Cette technique permettra de dégager des tendances et prédictions des couts des sinistres en se basant sur plusieurs variables explicatives. Nous allons commencer par identifier les variables prédictrices les plus pertinentes devant intervenir dans la modélisation du cout des sinistres. Par la suite nous allons modéliser par régression linéaire multiple les couts des sinistres. Et enfin nous présenterons une analyse de l’effet de la zone sur cout des sinistres dans cette compagnie d’assurances.
PARTIE I : TRAITEMENT DES DONNEES BRUTES
Cette partie présente le dictionnaire des données utilisées dans la présente analyse ainsi que les étapes de l’apurement de notre base de données. A cela, il faut ajouter que la méthodologie est basée sur une approche essentiellement descriptive en utilisant le logiciel R pour le traitement et l’analyse des données.
Dictionnaire des données de l’analyse
Le jeu de données s’appelera “Petrole”. Il comporte des variables tant quantitatives que qualitatives. Ces variables sont listées dans le tableau de dictionnaire de données ci-dessous :
NOM DES VARIABLES DESCRIPTION Cout Coût du sinistre Nocontrat Numero de contrat Exposition Durée du la signature du contrat Zone Commune ou quartier puissance Puissance du véhicule agevehicule Age du véhicule ageconducteur Age du conducteur bonus Coefficient de réduction-majoration marque Marque du véhicule carburant Type de carburant utilisé densite Densité de la population region La région nbre Le nombre de chauffeur du véhicule No Numéro du contrat d’assurance garantie La garantie donnée par le propriétaire du véhicule
my_tbl <- tibble::tribble(
~VARIABLE, ~NATURE, ~DESCRIPTION, ~MODALITES,
"Cout", "Quantitative", " Coût du sinistre", "Numérique entier",
"glucose", "Quantitative", "Concentration de glucose plasmatique ", " Numérique entier",
"pression_art", "Quantitative", " Pression artérielle diastolique (mm Hg)", "Numérique entier",
"triceps", "Quantitative", "Épaisseur du pli cutané du triceps (mm)", "Numérique entier",
"insuline", "Quantitative", " Insuline sérique(mu U/ml)", "Numérique entier",
"imc ", "Quantitative", "Indice de masse corporelle", "Numérique décimal",
"pedigree", "Quantitative", "Fonction pédigrée du diabète", "Numérique décimal",
"age", "Quantitative", " Âge (années)", "Numérique entier",
"diabete", "Qualitative", " Classe, apparition du diabète dans les cinq ans", "pos\nneg"
)
require(rhandsontable)
rhandsontable(my_tbl, rowHeaders = NULL,
digits = 3, useTypes = FALSE, search = FALSE,
width = NULL, height = NULL)Importation et transformation du jeu de donées
Lecture du jeu de données
Nous avons commencé par charger le jeu de données depuis son espace de stockage sur notre ordinateur avant de le lire. Cette base de données sera appelée « RegMul ». Par la suite, nous avons choisi d’en afficher un aperçu. Le jeu de données global contient 2765 observations reparties sur 15 variables dont 3 qualitatives et 12 quantitatives.
#Donner le chemin d'acces au jeu de données
setwd("D:/Insseds/Dataframe")
#Dataframe avec quanti & quali
RegPois <- read.table("actuarLog.csv",header=TRUE,sep=";",check.names=FALSE,stringsAsFactors = TRUE, na.strings = 0) # Lire le jeu de données
introduce(RegPois)## rows columns discrete_columns continuous_columns all_missing_columns
## 1 2765 15 3 12 0
## total_missing_values complete_rows total_observations memory_usage
## 1 516 2277 41475 193984summary(RegPois)## nocontrat exposition zone puissance agevehicule
## Min. : 217 Min. :0.008219 A:444 Min. : 4.000 Min. : 1.000
## 1st Qu.: 108786 1st Qu.:0.500000 B:300 1st Qu.: 5.000 1st Qu.: 3.000
## Median :1041406 Median :0.870000 C:760 Median : 6.000 Median : 6.000
## Mean : 778477 Mean :0.737445 D:619 Mean : 6.374 Mean : 6.624
## 3rd Qu.:1128643 3rd Qu.:1.000000 E:567 3rd Qu.: 7.000 3rd Qu.:10.000
## Max. :2054297 Max. :1.300000 F: 75 Max. :15.000 Max. :35.000
## NA's :173
## ageconducteur bonus marque carburant densite
## Min. :18.00 Min. : 50.00 Min. : 1.000 D:1434 Min. :11.00
## 1st Qu.:33.00 1st Qu.: 50.00 1st Qu.: 2.000 E:1331 1st Qu.:24.00
## Median :43.00 Median : 50.00 Median : 2.000 Median :52.00
## Mean :44.09 Mean : 61.21 Mean : 4.437 Mean :49.11
## 3rd Qu.:53.00 3rd Qu.: 68.00 3rd Qu.: 6.000 3rd Qu.:82.00
## Max. :99.00 Max. :165.00 Max. :14.000 Max. :94.00
##
## region nbre no garantie cout
## Min. :-1.00 Min. :1.000 Min. : 148 1RC:985 Min. : -3811.2
## 1st Qu.: 7.00 1st Qu.:1.000 1st Qu.: 28810 2DO:820 1st Qu.: 215.8
## Median :13.00 Median :1.000 Median : 38506 3VI:187 Median : 500.6
## Mean :10.22 Mean :1.712 Mean : 44531 4BG:761 Mean : 1190.8
## 3rd Qu.:13.00 3rd Qu.:2.000 3rd Qu.: 49055 5CO: 9 3rd Qu.: 1128.1
## Max. :13.00 Max. :7.000 Max. :104152 6CL: 3 Max. :152449.0
## NA's :61 NA's :282Recodage et transformation des variables
Nous avons dans cette partie procédé à un recodage des variables qualitatives dont les natures d’origine n’ont pas été lues.
RegPois$zone=factor(RegPois$zone,labels=c("Sud","Nord","Centre","Est","Cotière","Ouest"))
RegPois$marque=factor(RegPois$marque,labels=c("Toyota","Mazda","Ford","Renault","Nissan","Peugeot","Honda","Bmw","Range","Volkswagen","Suziki"))
RegPois$region=factor(RegPois$region,labels=c("Abidjan","Bouake","Yakro","korhogo","Odienne","Man","Daloa","Duekoue","San pedro","Lahou","Soubre", "Aboisso", "Bongouanou", "Daoukro"))
RegPois$carburant=factor(RegPois$carburant,labels=c("Diesel","Essence"))
attach(RegPois)
summary(nbre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 1.000 1.000 1.712 2.000 7.000RegPois$accident[nbre]## NULLRegPois$no=NULL
RegPois$nocontrat=NULL
summary(RegPois)## exposition zone puissance agevehicule
## Min. :0.008219 Sud :444 Min. : 4.000 Min. : 1.000
## 1st Qu.:0.500000 Nord :300 1st Qu.: 5.000 1st Qu.: 3.000
## Median :0.870000 Centre :760 Median : 6.000 Median : 6.000
## Mean :0.737445 Est :619 Mean : 6.374 Mean : 6.624
## 3rd Qu.:1.000000 Cotière:567 3rd Qu.: 7.000 3rd Qu.:10.000
## Max. :1.300000 Ouest : 75 Max. :15.000 Max. :35.000
## NA's :173
## ageconducteur bonus marque carburant densite
## Min. :18.00 Min. : 50.00 Mazda :755 Diesel :1434 Min. :11.00
## 1st Qu.:33.00 1st Qu.: 50.00 Toyota :683 Essence:1331 1st Qu.:24.00
## Median :43.00 Median : 50.00 Range :376 Median :52.00
## Mean :44.09 Mean : 61.21 Ford :289 Mean :49.11
## 3rd Qu.:53.00 3rd Qu.: 68.00 Nissan :188 3rd Qu.:82.00
## Max. :99.00 Max. :165.00 Peugeot:128 Max. :94.00
## (Other):346
## region nbre garantie cout
## Daoukro:1682 Min. :1.000 1RC:985 Min. : -3811.2
## Man : 114 1st Qu.:1.000 2DO:820 1st Qu.: 215.8
## Odienne: 102 Median :1.000 3VI:187 Median : 500.6
## Duekoue: 100 Mean :1.712 4BG:761 Mean : 1190.8
## Daloa : 88 3rd Qu.:2.000 5CO: 9 3rd Qu.: 1128.1
## (Other): 618 Max. :7.000 6CL: 3 Max. :152449.0
## NA's : 61 NA's :282Traitement des valeurs manquantes dans les données
Le résumé sommaire de ce jeu de données a permis de détecter qu’il existe des données manquantes matérialisées par des valeurs NAs. Ces valeurs doivent passer par une étape de traitement. Et cette section va servir à afficher les observations qui contiennent des valeurs manquantes. Elle srvira également à réaliser tout le traitement qui y convient.
Proportion des valeurs manquantes dans les données
Il existe 488 valeurs manquantes sur l’ensemble des données. Et le taux global de valeurs manquantes est de 17,6%.
nrow(RegPois[!complete.cases(RegPois),]) # Determiner le nombre de valeurs manquantes## [1] 488nrow(RegPois[!complete.cases(RegPois),])/nrow(RegPois)# déterminer le taux de valeurs manquantes## [1] 0.1764919Notre jeu de données utilise une memoire de 189,4Kb. En le visualisant, l’on se rend compte qu’en moyenne c’est 1,2% de données manquantes sur chacune des variables. la variable «cout» contient à elle seule 10% des données manquantes. Cette variable étant la variable à expliquer dans notre modelisation, nous allons procéder à un retraitement des données manquantes dans l’ensemble pour ne pas nous donner le luxe de les supprimer.
library(visdat)
vis_miss(RegPois)# Visulaiser le jeu de données entier avec les valeurs manquantes en pourcentage de NA pour chaque variable et global
vis_dat(RegPois)# Visualiser la position des données manquantes
library(naniar)
naniar::gg_miss_upset(RegPois)# savoir si les valeurs ne sont pas correllées
plot_intro(RegPois) 
plot_missing(RegPois)
Comparaison et sélection de méthodes d’imputation
Notre jeu de données brutes est hétérogène. Alors kNN est la méthode de complétion qui nous permettra d’imputer les données hétérogènes. Cette méthode est choisie en ce sens que le comportement des données est ensuite étudié lorsque la quantité de données manquantes augmente. L’erreur d’imputation des données quantitatives est calculée comme suit : l’on prendra la valeur absolue de la différence avec l’échantillon test. Pour les variables qualitatives, on utilisera la distance de Hamming.
# tests des méthodes d’imputation par des tests d'erreurs d'imputation pour retenir la meilleure
err.model<-function(RegPois,RegPois.model,ind.test){
err={}
for (i in sort(unique(ind.test[,2]))){
test=RegPois[ind.test[ind.test[,2]==i,1],i]
if (length(test)>0){
if (is.factor(RegPois[,i])){
#distance de Hamming (pourcentage de mauvais choix)
err=rbind(err,100*length(which(RegPois.model[,i]!=
RegPois[,i]))/length(test))
} else {
#moyenne de l’erreur en valeur absolue
err=rbind(err,mean(abs(as.numeric(
RegPois[ind.test[ind.test[,2]==i,1],i])-
as.numeric(RegPois.model[ind.test[
ind.test[,2]==i,1],i]))))
}}}
err
}Création de 10 à 80% de données manquantes artificielles
En se limitant toujours au cas MCAR, on va créer artificiellement de plus en plus de données manquantes à imputer : de 10 à 80%, ratio donné par TEST.RATIO. On initialisera la matrice d’erreurs, une ligne par ratio de données manquantes, une colonne par variable.
# Création de données manquantes de 10 à 80% de données manquantes crée artificiellement de plus en plus de données manquantes à imputer
tmp=1
TEST.RATIO=seq(0.1,0.8,by=0.1)
# initialisation des matrices d’erreur
err.kNN=matrix(NA,nrow=length(TEST.RATIO),ncol=13)
colnames(err.kNN)=names(RegPois)Imputation des valeurs manquantes
Nous avons d’abord créé l’échantillon test. Par la suite nous avons réalisé les imputations avec les 3 algoritmes. Le jeu de données de départ est resté le jeu d’apprentissage.
# Imputation de valeurs manquantes
for(test.ratio in TEST.RATIO){
# création de l’échantillon test
IND=which(!is.na(RegPois),arr.ind=TRUE)
ntest=ceiling(dim(RegPois)[1]*test.ratio)
ind.test=IND[sample(1:dim(IND)[1],ntest),]
RegPois.test=RegPois[ind.test]
RegPois.train=RegPois
RegPois.train[ind.test]=NA
# par kNN
RegPois.kNN=kNN(RegPois.train, k=50, imp_var=FALSE)
err.kNN[tmp,1:length(unique(ind.test[,2]))]=
err.model(RegPois,RegPois.kNN,ind.test)
}Traitement des valeurs abérrantes et extremes
Une seconde étape de l’exploration des données disponibles nous a conduit à traiter les valeurs aberrantes et extrêmes sur les variables quantitatives. Nous appelons valeur aberrante une valeur ou une observation qui est « distante » des autres observations effectuées sur le même phénomène, c’est-à-dire qu’elle contraste grandement avec les valeurs « normalement » mesurées. Leur présence dans les données peut conduire à des estimateurs de paramètres biaisés et, suite à la réalisation de tests statistiques, à une interprétation des résultats erronée. Pouvant être dû à plusieurs facteurs, nous avons pensé utile dans un premier temps de les detecter, puis de les imputer si elles existaient dans notre base de données. La détection desdites données est perceptible comme l’on peut le voir sur notre graphe. Nous avons utilisé les boite à moustache afin de visualiser les débordements observés.
Détection visuelle des valeurs aberrantes et extrêmes
library(rpart)
par(mfrow=c(2,4), mar=c(3,3,3,3))
boxplot(RegPois.kNN$exposition)
boxplot(RegPois.kNN$puissance)
boxplot(RegPois.kNN$agevehicule)
boxplot(RegPois.kNN$ageconducteur)
boxplot(RegPois.kNN$bonus)
boxplot(RegPois.kNN$densite)
boxplot(RegPois.kNN$cout)
Tous les graphiques présentant des points au-dessus ou en dessous de la boite sont des valeurs aberrantes. Seule la variable “glucose” ne présentent pas de valeurs aberrantes. Les autres variables qui présentent des données abérrantes seront retraitées.
Technique d’imputation de données abérrantes par winzorisation
Pour le traitement des données extrêmes nous avons utilisé la technique de Winzorisation en les ramenant dans les limites des bornes (inférieure et supérieure) des moutâches.
RegPois.kNN$puissance <- Winsorize(RegPois.kNN$puissance)
RegPois.kNN$agevehicule<- Winsorize(RegPois.kNN$agevehicule)
RegPois.kNN$ageconducteur<- Winsorize(RegPois.kNN$ageconducteur)
RegPois.kNN$bonus<- Winsorize(RegPois.kNN$bonus)
RegPois.kNN$cout<- Winsorize(RegPois.kNN$cout)
par(mfrow=c(2,3), mar=c(3,3,3,3))
boxplot(RegPois.kNN$puissance)# graph
boxplot(RegPois.kNN$agevehicule)
boxplot(RegPois.kNN$ageconducteur)
boxplot(RegPois.kNN$bonus)
boxplot(RegPois.kNN$cout)
par(mfrow=c(1,1), mar=c(0,0,0,0))
Toutes les valeurs abérrantes et extremes ont été traitées comme le montre les graphiques ci-dessus.
Statistique sommaire la variable cible
freq(RegPois.kNN$nbre)## n % val%
## 1 1473 53.3 53.3
## 2 895 32.4 32.4
## 3 248 9.0 9.0
## 4 111 4.0 4.0
## 5 18 0.7 0.7
## 6 13 0.5 0.5
## 7 7 0.3 0.3hist(RegPois.kNN$nbre)
plot(table(RegPois.kNN$nbre))
La distribution a une allure en « L », calculons maintenant la moyenne et la variance. En effet dans une distribution de poisson ces deux valeurs doivent être égales.
mean(RegPois.kNN$nbre)## [1] 1.688246var(RegPois.kNN$nbre)## [1] 0.8832373Ici la moyenne n’est pas très loin de la variance, on peut utiliser la régression de poisson. Si par contre la Variance est très supérieur à la moyennes: ce phénomène est appelé «sur-dispersion». Il peut se produire pour de nombreuses raisons, par exemple la nonindépendance des événements comptés. Pour faire face à la dispersion excessive, un modèle de régression quasi-Poisson peut être utilisé. Avec la fonction glm () et l’option quasipoisson.
Test d’adéquation à la loi de poisson
test d’adéquation avec une distribution de poisson
library(fitdistrplus)
#Estimation des paramètres par la méthode du maximum de vraisemblance
fpois <- fitdist(RegPois.kNN$nbre, "pois")
summary(fpois)## Fitting of the distribution ' pois ' by maximum likelihood
## Parameters :
## estimate Std. Error
## lambda 1.688246 0.02470985
## Loglikelihood: -3872.283 AIC: 7746.566 BIC: 7752.491#Tests d’adéquation
fpois <- fitdist(RegPois.kNN$nbre, "pois")
gofstat(fpois)## Chi-squared statistic: 153.343
## Degree of freedom of the Chi-squared distribution: 3
## Chi-squared p-value: 5.00684e-33
## Chi-squared table:
## obscounts theocounts
## <= 1 1473.00000 1373.94149
## <= 2 895.00000 728.35092
## <= 3 248.00000 409.87849
## <= 4 111.00000 172.99392
## > 4 38.00000 79.83517
##
## Goodness-of-fit criteria
## 1-mle-pois
## Akaike's Information Criterion 7746.566
## Bayesian Information Criterion 7752.491plot(fpois)
NB : p.value inférieur à 5%, les données ne suivent pas une loi de poisson
PARTIE II : MODELISATION DE LA PROBABILITE DE SURVENANCE DES ACCIDENTS
Jeux de données à utiliser pour la construction du modele
Nous allons utiliser le jeu de données avec les valeurs manquantes retraitées et débarrassé des valeurs abérrantes et extremes. RegLog.knn
Nous allons déterminer le jeu de données d’apprentissage sur les 2700 premières observations de notre jeu de données initial. RegLogTrain
Nous allons utiliser les 65 derbières observations de notre jeu de données de test. RegLogTest
df<-data.frame(RegPois.kNN)
# Echantillon d'apprentissage
RegPoisTrain<- df[-c(2701:2765), ]
# Echantillon de test, prediction
RegPoisTest<- df[-c(1:2700), ]Construction du modèle
Préselection des variables optimales du modele
Dans cette partie nous allons identifier les variables prédictrices les plus pertinentes devant intervenir dans la modélisation de la probabilité de survenance des accidents en fonction du profil de l’assuré.
L’algorithme de sélection automatique de variables à permis de detecter les variables exposition + zone + agevehicule + carburant + garantie + cout comme expliquant au mieux la survenue d’accidents.
library(MASS)
modePois <- stepAIC(glm(nbre~.,data=RegPoisTrain,family=poisson))## Start: AIC=7476.68
## nbre ~ exposition + zone + puissance + agevehicule + ageconducteur +
## bonus + marque + carburant + densite + region + garantie +
## cout
##
## Df Deviance AIC
## - region 13 1025.5 7461.6
## - zone 5 1018.6 7470.8
## - ageconducteur 1 1014.5 7474.7
## - bonus 1 1014.5 7474.7
## - densite 1 1014.6 7474.8
## - puissance 1 1014.7 7474.8
## <none> 1014.5 7476.7
## - cout 1 1018.9 7479.1
## - agevehicule 1 1020.5 7480.7
## - marque 10 1039.7 7481.8
## - exposition 1 1025.5 7485.7
## - carburant 1 1027.6 7487.7
## - garantie 5 1097.8 7549.9
##
## Step: AIC=7461.64
## nbre ~ exposition + zone + puissance + agevehicule + ageconducteur +
## bonus + marque + carburant + densite + garantie + cout
##
## Df Deviance AIC
## - zone 5 1029.6 7455.7
## - ageconducteur 1 1025.5 7459.7
## - densite 1 1025.6 7459.8
## - puissance 1 1025.7 7459.9
## <none> 1025.5 7461.6
## - bonus 1 1029.1 7463.3
## - cout 1 1029.3 7463.5
## - agevehicule 1 1032.2 7466.3
## - marque 10 1051.8 7467.9
## - exposition 1 1036.6 7470.8
## - carburant 1 1039.0 7473.2
## - garantie 5 1107.6 7533.7
##
## Step: AIC=7455.73
## nbre ~ exposition + puissance + agevehicule + ageconducteur +
## bonus + marque + carburant + densite + garantie + cout
##
## Df Deviance AIC
## - ageconducteur 1 1029.6 7453.7
## - densite 1 1029.6 7453.8
## - puissance 1 1029.8 7453.9
## <none> 1029.6 7455.7
## - cout 1 1033.4 7457.5
## - bonus 1 1033.9 7458.1
## - agevehicule 1 1036.2 7460.3
## - marque 10 1056.8 7463.0
## - exposition 1 1040.3 7464.5
## - carburant 1 1043.4 7467.5
## - garantie 5 1110.8 7527.0
##
## Step: AIC=7453.74
## nbre ~ exposition + puissance + agevehicule + bonus + marque +
## carburant + densite + garantie + cout
##
## Df Deviance AIC
## - densite 1 1029.6 7451.8
## - puissance 1 1029.8 7451.9
## <none> 1029.6 7453.7
## - cout 1 1033.4 7455.5
## - bonus 1 1035.3 7457.5
## - agevehicule 1 1036.2 7458.4
## - marque 10 1056.8 7461.0
## - exposition 1 1040.3 7462.5
## - carburant 1 1043.7 7465.8
## - garantie 5 1111.4 7525.6
##
## Step: AIC=7451.76
## nbre ~ exposition + puissance + agevehicule + bonus + marque +
## carburant + garantie + cout
##
## Df Deviance AIC
## - puissance 1 1029.8 7449.9
## <none> 1029.6 7451.8
## - cout 1 1033.4 7453.6
## - bonus 1 1035.3 7455.5
## - agevehicule 1 1036.2 7456.4
## - marque 10 1056.8 7459.0
## - exposition 1 1040.3 7460.5
## - carburant 1 1043.7 7463.8
## - garantie 5 1111.6 7523.7
##
## Step: AIC=7449.94
## nbre ~ exposition + agevehicule + bonus + marque + carburant +
## garantie + cout
##
## Df Deviance AIC
## <none> 1029.8 7449.9
## - cout 1 1033.6 7451.8
## - bonus 1 1035.4 7453.6
## - agevehicule 1 1036.3 7454.4
## - marque 10 1057.6 7457.8
## - exposition 1 1040.5 7458.7
## - carburant 1 1044.0 7462.2
## - garantie 5 1111.6 7521.7Conclusion : Si ∗, l’influence de Xj sur Y est significative, si ∗∗, elle est très significative et si ∗ ∗ ∗, elle est hautement significative.
Coefficients du modèle de regression
Dans le cadre d’un modèle logistique, généralement on ne présente pas les coefficients du modèle mais leur valeur exponentielle, cette dernière correspondant en effet à des odds ratio, également appelés rapports des cotes. Un odds ratio de 1 signifie l’absence d’effet. Un odds ratio largement supérieur à 1 correspond à une augmentation du phénomène étudié et un odds ratio largement inféieur à 1 correspond à une diminution du phénomène étudié.
L’analyse nous permet de détecter la présence d’effets. En effet, les variable exposition diminue la probabilité de survenance d’accident tandis que la variable agevehicule augmente la probabilité d’accidents.
La fonction tbl_regression de l’extension gtsummary, qui a recours en interne à broom.helpers, permet d’obtenir un tableau plus propre. Comme nous souhaitons afficher les odds ratios plutôt que les coefficients du modèle, on indiquera exponentiate = TRUE.
library(gtsummary)
tbl_regression(modePois, exponentiate = TRUE)| Characteristic | IRR1 | 95% CI1 | p-value |
|---|---|---|---|
| exposition | 1.19 | 1.07, 1.32 | 0.001 |
| agevehicule | 0.99 | 0.98, 1.00 | 0.011 |
| bonus | 1.00 | 1.00, 1.00 | 0.017 |
| marque | |||
| Toyota | — | — | |
| Mazda | 1.00 | 0.92, 1.09 | >0.9 |
| Ford | 1.10 | 0.99, 1.22 | 0.089 |
| Renault | 1.27 | 1.10, 1.46 | <0.001 |
| Nissan | 1.03 | 0.91, 1.17 | 0.6 |
| Peugeot | 1.05 | 0.90, 1.22 | 0.5 |
| Honda | 1.06 | 0.88, 1.27 | 0.5 |
| Bmw | 1.40 | 1.16, 1.68 | <0.001 |
| Range | 1.05 | 0.94, 1.17 | 0.4 |
| Volkswagen | 0.86 | 0.67, 1.09 | 0.2 |
| Suziki | 1.26 | 0.86, 1.77 | 0.2 |
| carburant | |||
| Diesel | — | — | |
| Essence | 0.89 | 0.84, 0.95 | <0.001 |
| garantie | |||
| 1RC | — | — | |
| 2DO | 1.07 | 0.99, 1.15 | 0.087 |
| 3VI | 1.03 | 0.91, 1.17 | 0.6 |
| 4BG | 0.77 | 0.71, 0.83 | <0.001 |
| 5CO | 1.98 | 1.36, 2.77 | <0.001 |
| 6CL | 0.84 | 0.26, 1.96 | 0.7 |
| cout | 1.00 | 1.00, 1.00 | 0.049 |
| 1 IRR = Incidence Rate Ratio, CI = Confidence Interval | |||
Identifier les variables ayant un effet significatif
Les p-values associées aux odds ratios nous indique si un odd ratio est significativement différent de 1, par rapport à la modalité de référence. Mais cela n’indique pas si globalement une variable a un effet significatif sur le modèle. Pour tester l’effet global sur un modèle, on peut avoir recours à la fonction drop1. Cette dernière va tour à tour supprimer chaque variable du modèle et réaliser une analyse de variance (ANOVA, voir fonction anova) pour voir si la variance change significativement.
un test d’anova permet de voir les variables les plus liées à la cible (ici le nombre d’accidents)
drop1(modePois, test = "Chisq")## Single term deletions
##
## Model:
## nbre ~ exposition + agevehicule + bonus + marque + carburant +
## garantie + cout
## Df Deviance AIC LRT Pr(>Chi)
## <none> 1029.8 7449.9
## exposition 1 1040.5 7458.7 10.724 0.0010578 **
## agevehicule 1 1036.3 7454.4 6.508 0.0107391 *
## bonus 1 1035.4 7453.6 5.631 0.0176420 *
## marque 10 1057.6 7457.8 27.857 0.0019026 **
## carburant 1 1044.0 7462.2 14.231 0.0001617 ***
## garantie 5 1111.6 7521.7 81.797 3.53e-16 ***
## cout 1 1033.6 7451.8 3.811 0.0509162 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modePois %>%
tbl_regression(exponentiate = TRUE) %>%
add_global_p(type = "II")| Characteristic | IRR1 | 95% CI1 | p-value |
|---|---|---|---|
| exposition | 1.19 | 1.07, 1.32 | 0.001 |
| agevehicule | 0.99 | 0.98, 1.00 | 0.011 |
| bonus | 1.00 | 1.00, 1.00 | 0.018 |
| marque | 0.002 | ||
| Toyota | — | — | |
| Mazda | 1.00 | 0.92, 1.09 | |
| Ford | 1.10 | 0.99, 1.22 | |
| Renault | 1.27 | 1.10, 1.46 | |
| Nissan | 1.03 | 0.91, 1.17 | |
| Peugeot | 1.05 | 0.90, 1.22 | |
| Honda | 1.06 | 0.88, 1.27 | |
| Bmw | 1.40 | 1.16, 1.68 | |
| Range | 1.05 | 0.94, 1.17 | |
| Volkswagen | 0.86 | 0.67, 1.09 | |
| Suziki | 1.26 | 0.86, 1.77 | |
| carburant | <0.001 | ||
| Diesel | — | — | |
| Essence | 0.89 | 0.84, 0.95 | |
| garantie | <0.001 | ||
| 1RC | — | — | |
| 2DO | 1.07 | 0.99, 1.15 | |
| 3VI | 1.03 | 0.91, 1.17 | |
| 4BG | 0.77 | 0.71, 0.83 | |
| 5CO | 1.98 | 1.36, 2.77 | |
| 6CL | 0.84 | 0.26, 1.96 | |
| cout | 1.00 | 1.00, 1.00 | 0.051 |
| 1 IRR = Incidence Rate Ratio, CI = Confidence Interval | |||
Par défaut, Anova utilise le type II et add_global_p le type III. Dans les deux cas, il est possible de préciser le type de test avec type = “II” ou type = “III”.
Dans le cas de notre exemple, un modèle simple sans interaction, le type de test ne change pas les résultats.
anova(modePois, type = "II",test="Chisq" )## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: nbre
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 2699 1194.4
## exposition 1 2.288 2698 1192.2 0.130418
## agevehicule 1 26.843 2697 1165.3 2.207e-07 ***
## bonus 1 8.436 2696 1156.9 0.003679 **
## marque 10 27.735 2686 1129.1 0.001991 **
## carburant 1 8.369 2685 1120.8 0.003816 **
## garantie 5 87.180 2680 1033.6 < 2.2e-16 ***
## cout 1 3.811 2679 1029.8 0.050916 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Représentation graphique du modèle
L’extension forestmodel propose de son côté une fonction forest_model qui, à partir d’un modèle, propose une représentation visuelle et tabulaire des coefficients.
library(forestmodel)
forest_model(modePois)
Représentation graphique des effets
L’extension effects propose une représentation graphique résumant les effets de chaque variable du modèle. Pour cela, il suffit d’appliquer la méthode plot au résultat de la fonction allEffects. Nous obtenons alors la figure ci-dessous.
library(effects)
library(ggeffects)
ggeffect(modePois, "bonus")## # Predicted counts of nbre
##
## bonus | Predicted | 95% CI
## --------------------------------
## 50 | 1.62 | [1.56, 1.68]
## 56 | 1.64 | [1.59, 1.70]
## 62 | 1.67 | [1.62, 1.72]
## 68 | 1.69 | [1.64, 1.75]
## 72 | 1.71 | [1.65, 1.78]
## 78 | 1.74 | [1.66, 1.82]
## 84 | 1.76 | [1.67, 1.86]
## 96 | 1.81 | [1.68, 1.96]plot(ggeffect(modePois, "bonus"))
ggeffect(modePois, "agevehicule")## # Predicted counts of nbre
##
## agevehicule | Predicted | 95% CI
## --------------------------------------
## 1 | 1.76 | [1.67, 1.85]
## 3 | 1.72 | [1.66, 1.79]
## 4 | 1.70 | [1.65, 1.76]
## 6 | 1.67 | [1.62, 1.72]
## 8 | 1.63 | [1.58, 1.69]
## 9 | 1.61 | [1.55, 1.68]
## 11 | 1.58 | [1.50, 1.66]
## 14 | 1.53 | [1.42, 1.64]plot(ggeffect(modePois, "agevehicule"))
ggeffect(modePois, "garantie")## # Predicted counts of nbre
##
## garantie | Predicted | 95% CI
## -----------------------------------
## 1RC | 1.75 | [1.67, 1.83]
## 2DO | 1.86 | [1.77, 1.97]
## 3VI | 1.81 | [1.62, 2.02]
## 4BG | 1.34 | [1.26, 1.43]
## 5CO | 3.46 | [2.44, 4.92]
## 6CL | 1.47 | [0.55, 3.92]plot(ggeffect(modePois, "garantie"))
ggeffect(modePois, "marque")## # Predicted counts of nbre
##
## marque | Predicted | 95% CI
## ----------------------------------
## Toyota | 1.60 | [1.51, 1.70]
## Mazda | 1.60 | [1.51, 1.69]
## Renault | 2.04 | [1.80, 2.31]
## Nissan | 1.65 | [1.47, 1.85]
## Peugeot | 1.68 | [1.46, 1.93]
## Honda | 1.70 | [1.43, 2.01]
## Bmw | 2.24 | [1.88, 2.68]
## Suziki | 2.02 | [1.42, 2.87]plot(ggeffect(modePois, "marque"))
ggeffect(modePois, "garantie")## # Predicted counts of nbre
##
## garantie | Predicted | 95% CI
## -----------------------------------
## 1RC | 1.75 | [1.67, 1.83]
## 2DO | 1.86 | [1.77, 1.97]
## 3VI | 1.81 | [1.62, 2.02]
## 4BG | 1.34 | [1.26, 1.43]
## 5CO | 3.46 | [2.44, 4.92]
## 6CL | 1.47 | [0.55, 3.92]plot(ggeffect(modePois, "garantie"))
Tableaux all-in-one
Nous avons déjà vu précédemment la fonction tbl_regression de l’extension gtsummary. Une vignette dédiée de l’extension explicite les possibilités de personnalisation des résultats. Par exemple add_global_p permet d’ajouter des p-valeurs globales pour chaque variable.
library(gtsummary)
theme_gtsummary_language("fr", decimal.mark = ",", big.mark = " ")
tbl_regression(modePois, exponentiate = TRUE) %>%
add_global_p(keep = TRUE)| Caractéristique | IRR1 | 95% IC1 | p-valeur |
|---|---|---|---|
| exposition | 1,19 | 1,07 – 1,32 | 0,001 |
| agevehicule | 0,99 | 0,98 – 1,00 | 0,011 |
| bonus | 1,00 | 1,00 – 1,00 | 0,018 |
| marque | 0,002 | ||
| Toyota | — | — | |
| Mazda | 1,00 | 0,92 – 1,09 | >0,9 |
| Ford | 1,10 | 0,99 – 1,22 | 0,089 |
| Renault | 1,27 | 1,10 – 1,46 | <0,001 |
| Nissan | 1,03 | 0,91 – 1,17 | 0,6 |
| Peugeot | 1,05 | 0,90 – 1,22 | 0,5 |
| Honda | 1,06 | 0,88 – 1,27 | 0,5 |
| Bmw | 1,40 | 1,16 – 1,68 | <0,001 |
| Range | 1,05 | 0,94 – 1,17 | 0,4 |
| Volkswagen | 0,86 | 0,67 – 1,09 | 0,2 |
| Suziki | 1,26 | 0,86 – 1,77 | 0,2 |
| carburant | <0,001 | ||
| Diesel | — | — | |
| Essence | 0,89 | 0,84 – 0,95 | <0,001 |
| garantie | <0,001 | ||
| 1RC | — | — | |
| 2DO | 1,07 | 0,99 – 1,15 | 0,087 |
| 3VI | 1,03 | 0,91 – 1,17 | 0,6 |
| 4BG | 0,77 | 0,71 – 0,83 | <0,001 |
| 5CO | 1,98 | 1,36 – 2,77 | <0,001 |
| 6CL | 0,84 | 0,26 – 1,96 | 0,7 |
| cout | 1,00 | 1,00 – 1,00 | 0,051 |
| 1 IRR = rapport de taux d'incidence, IC = intervalle de confiance | |||
Prédiction du nombre d’accident par assuré
Une manière de tester la qualité d’un modèle est le calcul d’une matrice de confusion, c’est-à-dire le tableau croisé des valeurs observées et celles des valeurs prédites en appliquant le modèle aux données d’origine.
La méthode predict avec l’argument type=“response” permet d’appliquer notre modèle logistique à un tableau de données et renvoie pour chaque individu le nombre probable d’accident qu’il fera.
nbre.pred <- predict(modePois, type = "response", newdata = RegPoisTest)
nbre.pred## 2701 2702 2703 2704 2705 2706 2707 2708
## 1.976726 1.797419 1.926553 1.539697 1.542811 2.197661 1.975066 1.635357
## 2709 2710 2711 2712 2713 2714 2715 2716
## 1.675777 1.759105 1.661569 1.776989 1.669527 1.608947 1.666333 1.256988
## 2717 2718 2719 2720 2721 2722 2723 2724
## 1.134439 1.575834 1.369126 1.151385 2.076309 1.896828 1.661282 1.947794
## 2725 2726 2727 2728 2729 2730 2731 2732
## 1.595399 1.715300 1.482622 1.508531 1.465761 1.724202 1.626546 1.252905
## 2733 2734 2735 2736 2737 2738 2739 2740
## 2.864903 1.179782 1.514496 1.874936 1.680503 1.520822 2.000141 1.606361
## 2741 2742 2743 2744 2745 2746 2747 2748
## 1.680338 1.457420 1.557749 1.724919 1.325949 1.730062 1.783971 1.558396
## 2749 2750 2751 2752 2753 2754 2755 2756
## 1.602227 1.938018 1.814596 1.379610 1.389421 1.169847 1.149510 1.391655
## 2757 2758 2759 2760 2761 2762 2763 2764
## 1.780540 1.939018 1.642181 2.014621 1.051307 1.394371 0.966604 1.471508
## 2765
## 1.399669L’extension finalfit fournit une fonction finalfit du type all-in-one qui calcule un tableau avec les tris croisés, les odds ratios univariés et un modèle multivarié.
Il faut d’abord définir la variable dépendante et les variables explicatives.
dep <- "nbre"
vars <- c("exposition", "agevehicule", "bonus","marque","carburant","garantie")Une première fonction summary_factorlist fournit un tableau descriptif avec, si l’option p = TRUE est indiquée, des tests de comparaisons (ici des tests du Chi²).
library(finalfit)
tab <- summary_factorlist(RegPoisTest, dep, vars, p = TRUE, add_dependent_label = TRUE)
knitr::kable(tab, row.names = FALSE)| Dependent: nbre | unit | value | p | |
|---|---|---|---|---|
| exposition | [0.0,1.0] | Mean (sd) | 1.5 (0.8) | 0.719 |
| agevehicule | [1.0,14.0] | Mean (sd) | 1.5 (0.8) | 0.895 |
| bonus | [50.0,95.0] | Mean (sd) | 1.5 (0.8) | 0.903 |
| marque | Toyota | Mean (sd) | 1.8 (1.1) | 0.186 |
| Mazda | Mean (sd) | 1.2 (0.4) | ||
| Ford | Mean (sd) | 1.0 (0.0) | ||
| Renault | Mean (sd) | 1.3 (0.6) | ||
| Peugeot | Mean (sd) | 1.0 (0.0) | ||
| Range | Mean (sd) | 1.9 (0.7) | ||
| carburant | Diesel | Mean (sd) | 1.6 (0.9) | 0.143 |
| Essence | Mean (sd) | 1.3 (0.5) | ||
| garantie | 1RC | Mean (sd) | 1.7 (0.9) | 0.123 |
| 2DO | Mean (sd) | 1.4 (0.6) | ||
| 3VI | Mean (sd) | 2.2 (1.3) | ||
| 4BG | Mean (sd) | 1.3 (0.7) |
Analyse des résidus
)
res.m<-rstudent(modePois)
plot(res.m,pch=15,cex=.5,ylab="Residus",ylim=c(-3,3))
abline(h=c(-2,0,2),lty=c(2,1,2))
En théorie 95% des résidus studentisés se trouvent dans l’intervalle [-2;2]. Ici on a visuellement beaucoup de résidus qui se trouvent dans cet intervalle. Ce qui est acceptable.
res.m<-rstudent(modePois)
sum(as.numeric(abs(res.m)<=3))/nrow(RegPoisTrain)*100 ## [1] 100