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")
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.
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.
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)
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
RegLog <- read.table("actuarLog.csv",header=TRUE,sep=";",check.names=FALSE,stringsAsFactors = TRUE, na.strings = 0) # Lire le jeu de données
introduce(RegLog)
## 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 193984
summary(RegLog)
## 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 :282
Nous avons dans cette partie procédé à un recodage des variables qualitatives dont les natures d’origine n’ont pas été lues.
RegLog$zone=factor(RegLog$zone,labels=c("Sud","Nord","Centre","Est","Cotière","Ouest"))
RegLog$marque=factor(RegLog$marque,labels=c("Toyota","Mazda","Ford","Renault","Nissan","Peugeot","Honda","Bmw","Range","Volkswagen","Suziki"))
RegLog$region=factor(RegLog$region,labels=c("Abidjan","Bouake","Yakro","korhogo","Odienne","Man","Daloa","Duekoue","San pedro","Lahou","Soubre", "Aboisso", "Bongouanou", "Daoukro"))
RegLog$carburant=factor(RegLog$carburant,labels=c("Diesel","Essence"))
attach(RegLog)
summary(nbre)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 1.000 1.000 1.712 2.000 7.000
RegLog$accident[nbre<=1]<-"Oui"
RegLog$accident[nbre>1]<-"non"
RegLog$accident=factor(RegLog$accident)
RegLog$no=NULL
RegLog$nbre=NULL
RegLog$nocontrat=NULL
summary(RegLog)
## 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 garantie cout accident
## Daoukro:1682 1RC:985 Min. : -3811.2 non:1314
## Man : 114 2DO:820 1st Qu.: 215.8 Oui:1451
## Odienne: 102 3VI:187 Median : 500.6
## Duekoue: 100 4BG:761 Mean : 1190.8
## Daloa : 88 5CO: 9 3rd Qu.: 1128.1
## (Other): 618 6CL: 3 Max. :152449.0
## NA's : 61 NA's :282
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.
Il existe 488 valeurs manquantes sur l’ensemble des données. Et le taux global de valeurs manquantes est de 17,6%.
nrow(RegLog[!complete.cases(RegLog),]) # Determiner le nombre de valeurs manquantes
## [1] 488
nrow(RegLog[!complete.cases(RegLog),])/nrow(RegLog)# déterminer le taux de valeurs manquantes
## [1] 0.1764919
Notre 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(RegLog)# Visulaiser le jeu de données entier avec les valeurs manquantes en pourcentage de NA pour chaque variable et global
vis_dat(RegLog)# Visualiser la position des données manquantes
library(naniar)
naniar::gg_miss_upset(RegLog)# savoir si les valeurs ne sont pas correllées
plot_intro(RegLog)
plot_missing(RegLog)
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(RegLog,RegLog.model,ind.test){
err={}
for (i in sort(unique(ind.test[,2]))){
test=RegLog[ind.test[ind.test[,2]==i,1],i]
if (length(test)>0){
if (is.factor(RegLog[,i])){
#distance de Hamming (pourcentage de mauvais choix)
err=rbind(err,100*length(which(RegLog.model[,i]!=
RegLog[,i]))/length(test))
} else {
#moyenne de l’erreur en valeur absolue
err=rbind(err,mean(abs(as.numeric(
RegLog[ind.test[ind.test[,2]==i,1],i])-
as.numeric(RegLog.model[ind.test[
ind.test[,2]==i,1],i]))))
}}}
err
}
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(RegLog)
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(RegLog),arr.ind=TRUE)
ntest=ceiling(dim(RegLog)[1]*test.ratio)
ind.test=IND[sample(1:dim(IND)[1],ntest),]
RegLog.test=RegLog[ind.test]
RegLog.train=RegLog
RegLog.train[ind.test]=NA
# par kNN
RegLog.kNN=kNN(RegLog.train, k=50, imp_var=FALSE)
err.kNN[tmp,1:length(unique(ind.test[,2]))]=
err.model(RegLog,RegLog.kNN,ind.test)
}
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.
library(rpart)
par(mfrow=c(2,4), mar=c(3,3,3,3))
boxplot(RegLog.kNN$exposition)
boxplot(RegLog.kNN$puissance)
boxplot(RegLog.kNN$agevehicule)
boxplot(RegLog.kNN$ageconducteur)
boxplot(RegLog.kNN$bonus)
boxplot(RegLog.kNN$densite)
boxplot(RegLog.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.
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.
RegLog.kNN$puissance <- Winsorize(RegLog.kNN$puissance)
RegLog.kNN$agevehicule<- Winsorize(RegLog.kNN$agevehicule)
RegLog.kNN$ageconducteur<- Winsorize(RegLog.kNN$ageconducteur)
RegLog.kNN$bonus<- Winsorize(RegLog.kNN$bonus)
RegLog.kNN$cout<- Winsorize(RegLog.kNN$cout)
par(mfrow=c(2,3), mar=c(3,3,3,3))
boxplot(RegLog.kNN$puissance)# graph
boxplot(RegLog.kNN$agevehicule)
boxplot(RegLog.kNN$ageconducteur)
boxplot(RegLog.kNN$bonus)
boxplot(RegLog.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.
freq(RegLog$accident)
## n % val%
## non 1314 47.5 47.5
## Oui 1451 52.5 52.5
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(RegLog.kNN)
# Echantillon d'apprentissage
RegLogTrain<- df[-c(2701:2765), ]
# Echantillon de test, prediction
RegLogTest<- df[-c(1:2700), ]
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)
modeLog <- stepAIC(glm(accident~.,data=RegLogTrain,family=binomial))
## Start: AIC=3362.65
## accident ~ exposition + zone + puissance + agevehicule + ageconducteur +
## bonus + marque + carburant + densite + region + garantie +
## cout
##
## Df Deviance AIC
## - zone 5 3286.4 3360.4
## - ageconducteur 1 3279.0 3361.0
## - densite 1 3279.1 3361.1
## - puissance 1 3279.5 3361.5
## <none> 3278.7 3362.7
## - bonus 1 3282.2 3364.2
## - region 13 3307.1 3365.1
## - cout 1 3290.5 3372.5
## - marque 10 3308.5 3372.5
## - agevehicule 1 3293.1 3375.1
## - carburant 1 3295.1 3377.1
## - exposition 1 3296.0 3378.0
## - garantie 5 3579.1 3653.1
##
## Step: AIC=3360.35
## accident ~ exposition + puissance + agevehicule + ageconducteur +
## bonus + marque + carburant + densite + region + garantie +
## cout
##
## Df Deviance AIC
## - ageconducteur 1 3286.7 3358.7
## - puissance 1 3287.4 3359.4
## - densite 1 3287.6 3359.6
## <none> 3286.4 3360.4
## - region 13 3314.1 3362.1
## - bonus 1 3290.1 3362.1
## - cout 1 3298.3 3370.3
## - marque 10 3316.9 3370.9
## - agevehicule 1 3300.7 3372.7
## - carburant 1 3303.1 3375.1
## - exposition 1 3303.2 3375.2
## - garantie 5 3583.7 3647.7
##
## Step: AIC=3358.71
## accident ~ exposition + puissance + agevehicule + bonus + marque +
## carburant + densite + region + garantie + cout
##
## Df Deviance AIC
## - puissance 1 3287.9 3357.9
## - densite 1 3288.0 3358.0
## <none> 3286.7 3358.7
## - bonus 1 3290.1 3360.1
## - region 13 3315.2 3361.2
## - cout 1 3298.9 3368.9
## - marque 10 3317.1 3369.1
## - agevehicule 1 3300.9 3370.9
## - carburant 1 3303.1 3373.1
## - exposition 1 3304.2 3374.2
## - garantie 5 3587.2 3649.2
##
## Step: AIC=3357.86
## accident ~ exposition + agevehicule + bonus + marque + carburant +
## densite + region + garantie + cout
##
## Df Deviance AIC
## - densite 1 3289.1 3357.1
## <none> 3287.9 3357.9
## - bonus 1 3291.2 3359.2
## - region 13 3316.9 3360.9
## - cout 1 3300.0 3368.0
## - marque 10 3318.3 3368.3
## - agevehicule 1 3301.5 3369.5
## - carburant 1 3304.7 3372.7
## - exposition 1 3305.2 3373.2
## - garantie 5 3587.3 3647.3
##
## Step: AIC=3357.1
## accident ~ exposition + agevehicule + bonus + marque + carburant +
## region + garantie + cout
##
## Df Deviance AIC
## <none> 3289.1 3357.1
## - bonus 1 3292.2 3358.2
## - region 13 3318.1 3360.1
## - marque 10 3319.2 3367.2
## - cout 1 3301.6 3367.6
## - agevehicule 1 3302.6 3368.6
## - carburant 1 3305.7 3371.7
## - exposition 1 3305.9 3371.9
## - garantie 5 3589.5 3647.5
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.
library(mfx)
logitmfx(accident~.,data=RegLogTrain)
## Call:
## logitmfx(formula = accident ~ ., data = RegLogTrain)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## exposition -1.5743e-01 3.7974e-02 -4.1459 3.385e-05 ***
## zoneNord -2.7018e-02 4.3519e-02 -0.6208 0.5347078
## zoneCentre -5.8307e-02 3.4301e-02 -1.6998 0.0891600 .
## zoneEst -8.4192e-02 3.6688e-02 -2.2948 0.0217432 *
## zoneCotière -1.7266e-02 3.7808e-02 -0.4567 0.6479089
## zoneOuest -3.3768e-02 7.2031e-02 -0.4688 0.6392098
## puissance -6.2894e-03 7.0245e-03 -0.8954 0.3705998
## agevehicule 1.1022e-02 2.9106e-03 3.7869 0.0001525 ***
## ageconducteur -6.0775e-04 1.0020e-03 -0.6066 0.5441396
## bonus -2.5176e-03 1.3378e-03 -1.8818 0.0598587 .
## marqueMazda -1.3465e-02 2.9258e-02 -0.4602 0.6453480
## marqueFord -9.7362e-02 3.9874e-02 -2.4418 0.0146156 *
## marqueRenault -1.3635e-01 5.3333e-02 -2.5566 0.0105711 *
## marqueNissan -5.9250e-02 4.6636e-02 -1.2705 0.2039150
## marquePeugeot -7.9857e-02 5.3543e-02 -1.4915 0.1358396
## marqueHonda -2.4882e-02 6.6427e-02 -0.3746 0.7079741
## marqueBmw -2.1897e-01 6.8026e-02 -3.2189 0.0012869 **
## marqueRange 3.9316e-02 3.8499e-02 1.0212 0.3071524
## marqueVolkswagen 1.2007e-01 7.6097e-02 1.5778 0.1146157
## marqueSuziki -1.3364e-01 1.4453e-01 -0.9246 0.3551682
## carburantEssence 9.0590e-02 2.2303e-02 4.0618 4.870e-05 ***
## densite -2.9099e-04 4.2170e-04 -0.6900 0.4901783
## regionBouake 7.8899e-02 9.1593e-02 0.8614 0.3890126
## regionYakro -2.2915e-01 8.1559e-02 -2.8097 0.0049589 **
## regionkorhogo -1.7625e-01 8.6342e-02 -2.0413 0.0412170 *
## regionOdienne -1.2265e-01 8.6908e-02 -1.4113 0.1581633
## regionMan -5.5431e-02 8.8448e-02 -0.6267 0.5308499
## regionDaloa -1.9254e-01 8.7031e-02 -2.2123 0.0269470 *
## regionDuekoue -1.5986e-01 8.9288e-02 -1.7904 0.0733968 .
## regionSan pedro 4.5005e-02 1.0879e-01 0.4137 0.6791111
## regionLahou -9.9749e-02 1.0007e-01 -0.9968 0.3188529
## regionSoubre -9.9699e-02 1.0448e-01 -0.9543 0.3399449
## regionAboisso -1.9675e-01 1.0122e-01 -1.9437 0.0519262 .
## regionBongouanou 3.3307e-02 1.1426e-01 0.2915 0.7706804
## regionDaoukro -1.2834e-01 8.0990e-02 -1.5846 0.1130623
## garantie2DO -1.6944e-01 2.5775e-02 -6.5739 4.901e-11 ***
## garantie3VI -7.4576e-02 4.4464e-02 -1.6772 0.0934996 .
## garantie4BG 2.9493e-01 2.3664e-02 12.4631 < 2.2e-16 ***
## garantie5CO -4.3321e-01 1.0316e-01 -4.1992 2.679e-05 ***
## garantie6CL 7.8813e-02 2.9282e-01 0.2692 0.7878117
## cout -4.3225e-05 1.2643e-05 -3.4189 0.0006287 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "zoneNord" "zoneCentre" "zoneEst" "zoneCotière"
## [5] "zoneOuest" "marqueMazda" "marqueFord" "marqueRenault"
## [9] "marqueNissan" "marquePeugeot" "marqueHonda" "marqueBmw"
## [13] "marqueRange" "marqueVolkswagen" "marqueSuziki" "carburantEssence"
## [17] "regionBouake" "regionYakro" "regionkorhogo" "regionOdienne"
## [21] "regionMan" "regionDaloa" "regionDuekoue" "regionSan pedro"
## [25] "regionLahou" "regionSoubre" "regionAboisso" "regionBongouanou"
## [29] "regionDaoukro" "garantie2DO" "garantie3VI" "garantie4BG"
## [33] "garantie5CO" "garantie6CL"
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(modeLog, exponentiate = TRUE)
| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| exposition | 0.54 | 0.40, 0.73 | <0.001 |
| agevehicule | 1.04 | 1.02, 1.07 | <0.001 |
| bonus | 0.99 | 0.98, 1.00 | 0.080 |
| marque | |||
| Toyota | — | — | |
| Mazda | 0.95 | 0.76, 1.19 | 0.7 |
| Ford | 0.67 | 0.49, 0.91 | 0.012 |
| Renault | 0.58 | 0.38, 0.89 | 0.014 |
| Nissan | 0.78 | 0.54, 1.12 | 0.2 |
| Peugeot | 0.72 | 0.47, 1.10 | 0.13 |
| Honda | 0.88 | 0.53, 1.46 | 0.6 |
| Bmw | 0.37 | 0.20, 0.67 | 0.001 |
| Range | 1.11 | 0.82, 1.50 | 0.5 |
| Volkswagen | 1.58 | 0.83, 3.08 | 0.2 |
| Suziki | 0.54 | 0.17, 1.74 | 0.3 |
| carburant | |||
| Diesel | — | — | |
| Essence | 1.43 | 1.21, 1.70 | <0.001 |
| region | |||
| Abidjan | — | — | |
| Bouake | 1.36 | 0.64, 2.91 | 0.4 |
| Yakro | 0.40 | 0.19, 0.82 | 0.013 |
| korhogo | 0.48 | 0.23, 0.99 | 0.048 |
| Odienne | 0.62 | 0.31, 1.24 | 0.2 |
| Man | 0.80 | 0.40, 1.58 | 0.5 |
| Daloa | 0.44 | 0.21, 0.92 | 0.030 |
| Duekoue | 0.55 | 0.26, 1.14 | 0.11 |
| San pedro | 1.20 | 0.50, 2.89 | 0.7 |
| Lahou | 0.68 | 0.31, 1.49 | 0.3 |
| Soubre | 0.68 | 0.30, 1.56 | 0.4 |
| Aboisso | 0.42 | 0.18, 1.00 | 0.050 |
| Bongouanou | 1.15 | 0.47, 2.86 | 0.8 |
| Daoukro | 0.59 | 0.30, 1.13 | 0.11 |
| garantie | |||
| 1RC | — | — | |
| 2DO | 0.52 | 0.42, 0.63 | <0.001 |
| 3VI | 0.75 | 0.53, 1.06 | 0.11 |
| 4BG | 3.59 | 2.87, 4.50 | <0.001 |
| 5CO | 0.10 | 0.01, 0.59 | 0.035 |
| 6CL | 1.33 | 0.12, 28.9 | 0.8 |
| cout | 1.00 | 1.00, 1.00 | <0.001 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
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.
Si l’on a recours à tbl_regression, on peut facilement ajouter les p-valeurs globales avec add_global_p.
drop1(modeLog, test = "Chisq")
## Single term deletions
##
## Model:
## accident ~ exposition + agevehicule + bonus + marque + carburant +
## region + garantie + cout
## Df Deviance AIC LRT Pr(>Chi)
## <none> 3289.1 3357.1
## exposition 1 3305.9 3371.9 16.783 4.191e-05 ***
## agevehicule 1 3302.6 3368.6 13.476 0.0002416 ***
## bonus 1 3292.2 3358.2 3.079 0.0792984 .
## marque 10 3319.2 3367.2 30.074 0.0008331 ***
## carburant 1 3305.7 3371.7 16.654 4.486e-05 ***
## region 13 3318.1 3360.1 29.023 0.0064975 **
## garantie 5 3589.5 3647.5 300.452 < 2.2e-16 ***
## cout 1 3301.6 3367.6 12.484 0.0004105 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modeLog %>%
tbl_regression(exponentiate = TRUE) %>%
add_global_p(type = "II")
| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| exposition | 0.54 | 0.40, 0.73 | <0.001 |
| agevehicule | 1.04 | 1.02, 1.07 | <0.001 |
| bonus | 0.99 | 0.98, 1.00 | 0.079 |
| marque | <0.001 | ||
| Toyota | — | — | |
| Mazda | 0.95 | 0.76, 1.19 | |
| Ford | 0.67 | 0.49, 0.91 | |
| Renault | 0.58 | 0.38, 0.89 | |
| Nissan | 0.78 | 0.54, 1.12 | |
| Peugeot | 0.72 | 0.47, 1.10 | |
| Honda | 0.88 | 0.53, 1.46 | |
| Bmw | 0.37 | 0.20, 0.67 | |
| Range | 1.11 | 0.82, 1.50 | |
| Volkswagen | 1.58 | 0.83, 3.08 | |
| Suziki | 0.54 | 0.17, 1.74 | |
| carburant | <0.001 | ||
| Diesel | — | — | |
| Essence | 1.43 | 1.21, 1.70 | |
| region | 0.006 | ||
| Abidjan | — | — | |
| Bouake | 1.36 | 0.64, 2.91 | |
| Yakro | 0.40 | 0.19, 0.82 | |
| korhogo | 0.48 | 0.23, 0.99 | |
| Odienne | 0.62 | 0.31, 1.24 | |
| Man | 0.80 | 0.40, 1.58 | |
| Daloa | 0.44 | 0.21, 0.92 | |
| Duekoue | 0.55 | 0.26, 1.14 | |
| San pedro | 1.20 | 0.50, 2.89 | |
| Lahou | 0.68 | 0.31, 1.49 | |
| Soubre | 0.68 | 0.30, 1.56 | |
| Aboisso | 0.42 | 0.18, 1.00 | |
| Bongouanou | 1.15 | 0.47, 2.86 | |
| Daoukro | 0.59 | 0.30, 1.13 | |
| garantie | <0.001 | ||
| 1RC | — | — | |
| 2DO | 0.52 | 0.42, 0.63 | |
| 3VI | 0.75 | 0.53, 1.06 | |
| 4BG | 3.59 | 2.87, 4.50 | |
| 5CO | 0.10 | 0.01, 0.59 | |
| 6CL | 1.33 | 0.12, 28.9 | |
| cout | 1.00 | 1.00, 1.00 | <0.001 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
modeLog %>%
tbl_regression(exponentiate = TRUE) %>%
add_global_p(type = "II")
| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| exposition | 0.54 | 0.40, 0.73 | <0.001 |
| agevehicule | 1.04 | 1.02, 1.07 | <0.001 |
| bonus | 0.99 | 0.98, 1.00 | 0.079 |
| marque | <0.001 | ||
| Toyota | — | — | |
| Mazda | 0.95 | 0.76, 1.19 | |
| Ford | 0.67 | 0.49, 0.91 | |
| Renault | 0.58 | 0.38, 0.89 | |
| Nissan | 0.78 | 0.54, 1.12 | |
| Peugeot | 0.72 | 0.47, 1.10 | |
| Honda | 0.88 | 0.53, 1.46 | |
| Bmw | 0.37 | 0.20, 0.67 | |
| Range | 1.11 | 0.82, 1.50 | |
| Volkswagen | 1.58 | 0.83, 3.08 | |
| Suziki | 0.54 | 0.17, 1.74 | |
| carburant | <0.001 | ||
| Diesel | — | — | |
| Essence | 1.43 | 1.21, 1.70 | |
| region | 0.006 | ||
| Abidjan | — | — | |
| Bouake | 1.36 | 0.64, 2.91 | |
| Yakro | 0.40 | 0.19, 0.82 | |
| korhogo | 0.48 | 0.23, 0.99 | |
| Odienne | 0.62 | 0.31, 1.24 | |
| Man | 0.80 | 0.40, 1.58 | |
| Daloa | 0.44 | 0.21, 0.92 | |
| Duekoue | 0.55 | 0.26, 1.14 | |
| San pedro | 1.20 | 0.50, 2.89 | |
| Lahou | 0.68 | 0.31, 1.49 | |
| Soubre | 0.68 | 0.30, 1.56 | |
| Aboisso | 0.42 | 0.18, 1.00 | |
| Bongouanou | 1.15 | 0.47, 2.86 | |
| Daoukro | 0.59 | 0.30, 1.13 | |
| garantie | <0.001 | ||
| 1RC | — | — | |
| 2DO | 0.52 | 0.42, 0.63 | |
| 3VI | 0.75 | 0.53, 1.06 | |
| 4BG | 3.59 | 2.87, 4.50 | |
| 5CO | 0.10 | 0.01, 0.59 | |
| 6CL | 1.33 | 0.12, 28.9 | |
| cout | 1.00 | 1.00, 1.00 | <0.001 |
| 1 OR = Odds 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(modeLog, type = "II",test="Chisq" )
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: accident
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 2699 3734.0
## exposition 1 3.512 2698 3730.5 0.0609156 .
## agevehicule 1 53.519 2697 3676.9 2.561e-13 ***
## bonus 1 5.529 2696 3671.4 0.0187017 *
## marque 10 21.316 2686 3650.1 0.0189978 *
## carburant 1 6.108 2685 3644.0 0.0134536 *
## region 13 28.188 2672 3615.8 0.0085180 **
## garantie 5 314.225 2667 3301.6 < 2.2e-16 ***
## cout 1 12.484 2666 3289.1 0.0004105 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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(modeLog)
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(modeLog, "carburant")
## # Predicted probabilities of accident
##
## carburant | Predicted | 95% CI
## ------------------------------------
## Diesel | 0.50 | [0.47, 0.52]
## Essence | 0.58 | [0.55, 0.61]
plot(ggeffect(modeLog, "carburant"))
ggeffect(modeLog, "agevehicule")
## # Predicted probabilities of accident
##
## agevehicule | Predicted | 95% CI
## --------------------------------------
## 1 | 0.48 | [0.45, 0.52]
## 3 | 0.50 | [0.47, 0.53]
## 4 | 0.51 | [0.49, 0.54]
## 6 | 0.53 | [0.51, 0.55]
## 8 | 0.56 | [0.53, 0.58]
## 10 | 0.58 | [0.55, 0.61]
## 12 | 0.60 | [0.56, 0.63]
## 15 | 0.63 | [0.58, 0.68]
plot(ggeffect(modeLog, "agevehicule"))
ggeffect(modeLog, "garantie")
## # Predicted probabilities of accident
##
## garantie | Predicted | 95% CI
## -----------------------------------
## 1RC | 0.51 | [0.47, 0.54]
## 2DO | 0.35 | [0.31, 0.38]
## 3VI | 0.43 | [0.36, 0.51]
## 4BG | 0.79 | [0.75, 0.81]
## 5CO | 0.09 | [0.01, 0.46]
## 6CL | 0.58 | [0.11, 0.94]
plot(ggeffect(modeLog, "garantie"))
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 la probabilité qu’il ait vécu le phénomène étudié c’est à dire la probabilité qu’il ait fait un accident ou non.
accident.pred <- predict(modeLog, type = "response", newdata = RegLogTest)
accident.pred
## 2701 2702 2703 2704 2705 2706 2707 2708
## 0.2807666 0.6210446 0.3884624 0.7023904 0.7060361 0.3257115 0.4303143 0.6247092
## 2709 2710 2711 2712 2713 2714 2715 2716
## 0.6044394 0.4766454 0.4803311 0.4394789 0.4037112 0.5399281 0.4148607 0.8206175
## 2717 2718 2719 2720 2721 2722 2723 2724
## 0.6440771 0.4086033 0.7730816 0.8489088 0.3313475 0.4961120 0.4873839 0.3819705
## 2725 2726 2727 2728 2729 2730 2731 2732
## 0.5598447 0.5255412 0.6399736 0.5977221 0.6033096 0.3826614 0.4784740 0.7878574
## 2733 2734 2735 2736 2737 2738 2739 2740
## 0.1206017 0.8468422 0.6785346 0.3180004 0.4074175 0.6180281 0.2568936 0.4558668
## 2741 2742 2743 2744 2745 2746 2747 2748
## 0.4078422 0.6499106 0.5495729 0.3811674 0.8187420 0.5536994 0.4591944 0.5741604
## 2749 2750 2751 2752 2753 2754 2755 2756
## 0.4203554 0.2850930 0.2258151 0.7724285 0.6592625 0.8405080 0.8351223 0.8263091
## 2757 2758 2759 2760 2761 2762 2763 2764
## 0.4536126 0.3407799 0.7186640 0.5361674 0.8559335 0.4706632 0.9180297 0.7808873
## 2765
## 0.7577800
Puisque notre variable étudiée est de type binaire. Nous devons donc transformer nos probabilités prédites en une variable du type « oui / non ». Usuellement, les probabilités prédites seront réunies en deux groupes selon qu’elles soient supérieures ou inférieures à la moitié. On peut ainsi dresser le tableau de contingence entre les labels prédits et les vrais labels.
table(accident.pred > 0.5, RegLogTest$accident)
##
## non Oui
## FALSE 12 18
## TRUE 15 20
Nous avons donc 29 (15+14) prédictions incorrectes sur un total de 65, soit un taux de mauvais classement de 44,6%.
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(modeLog, exponentiate = TRUE) %>%
add_global_p(keep = TRUE)
| Caractéristique | OR1 | 95% IC1 | p-valeur |
|---|---|---|---|
| exposition | 0,54 | 0,40 – 0,73 | <0,001 |
| agevehicule | 1,04 | 1,02 – 1,07 | <0,001 |
| bonus | 0,99 | 0,98 – 1,00 | 0,079 |
| marque | <0,001 | ||
| Toyota | — | — | |
| Mazda | 0,95 | 0,76 – 1,19 | 0,7 |
| Ford | 0,67 | 0,49 – 0,91 | 0,012 |
| Renault | 0,58 | 0,38 – 0,89 | 0,014 |
| Nissan | 0,78 | 0,54 – 1,12 | 0,2 |
| Peugeot | 0,72 | 0,47 – 1,10 | 0,13 |
| Honda | 0,88 | 0,53 – 1,46 | 0,6 |
| Bmw | 0,37 | 0,20 – 0,67 | 0,001 |
| Range | 1,11 | 0,82 – 1,50 | 0,5 |
| Volkswagen | 1,58 | 0,83 – 3,08 | 0,2 |
| Suziki | 0,54 | 0,17 – 1,74 | 0,3 |
| carburant | <0,001 | ||
| Diesel | — | — | |
| Essence | 1,43 | 1,21 – 1,70 | <0,001 |
| region | 0,006 | ||
| Abidjan | — | — | |
| Bouake | 1,36 | 0,64 – 2,91 | 0,4 |
| Yakro | 0,40 | 0,19 – 0,82 | 0,013 |
| korhogo | 0,48 | 0,23 – 0,99 | 0,048 |
| Odienne | 0,62 | 0,31 – 1,24 | 0,2 |
| Man | 0,80 | 0,40 – 1,58 | 0,5 |
| Daloa | 0,44 | 0,21 – 0,92 | 0,030 |
| Duekoue | 0,55 | 0,26 – 1,14 | 0,11 |
| San pedro | 1,20 | 0,50 – 2,89 | 0,7 |
| Lahou | 0,68 | 0,31 – 1,49 | 0,3 |
| Soubre | 0,68 | 0,30 – 1,56 | 0,4 |
| Aboisso | 0,42 | 0,18 – 1,00 | 0,050 |
| Bongouanou | 1,15 | 0,47 – 2,86 | 0,8 |
| Daoukro | 0,59 | 0,30 – 1,13 | 0,11 |
| garantie | <0,001 | ||
| 1RC | — | — | |
| 2DO | 0,52 | 0,42 – 0,63 | <0,001 |
| 3VI | 0,75 | 0,53 – 1,06 | 0,11 |
| 4BG | 3,59 | 2,87 – 4,50 | <0,001 |
| 5CO | 0,10 | 0,01 – 0,59 | 0,035 |
| 6CL | 1,33 | 0,12 – 28,9 | 0,8 |
| cout | 1,00 | 1,00 – 1,00 | <0,001 |
| 1 OR = rapport de cotes, IC = intervalle de confiance | |||
RegLogTrain %>%
dplyr::select(accident, zone, agevehicule, cout) %>%
tbl_summary(by = "accident", percent = "row") %>%
add_overall(last = TRUE) %>%
add_p()
| Caractéristique | non, N = 1 2721 | Oui, N = 1 4281 | Total, N = 2 7001 | p-valeur2 |
|---|---|---|---|---|
| zone | 0,4 | |||
| Sud | 176 (43%) | 238 (57%) | 414 (100%) | |
| Nord | 122 (46%) | 146 (54%) | 268 (100%) | |
| Centre | 382 (48%) | 409 (52%) | 791 (100%) | |
| Est | 294 (49%) | 302 (51%) | 596 (100%) | |
| Cotière | 264 (47%) | 297 (53%) | 561 (100%) | |
| Ouest | 34 (49%) | 36 (51%) | 70 (100%) | |
| agevehicule | 5,0 (3,0 – 7,0) | 6,0 (3,0 – 10,0) | 6,0 (3,0 – 9,0) | <0,001 |
| cout | 667 (286 – 1 128) | 423 (234 – 1 056) | 516 (251 – 1 128) | <0,001 |
| 1 n (%); Médiane (EI) | ||||
| 2 test du khi-deux d'indépendance; test de Wilcoxon-Mann-Whitney | ||||
L’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 <- "accident"
vars <- c("exposition", "puissance", "cout", "bonus","marque","region","zone")
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(RegLogTest, dep, vars, p = TRUE, add_dependent_label = TRUE)
knitr::kable(tab, row.names = FALSE)
| Dependent: accident | non | Oui | p | |
|---|---|---|---|---|
| exposition | Mean (SD) | 0.7 (0.4) | 0.7 (0.4) | 0.545 |
| puissance | Mean (SD) | 6.9 (2.0) | 6.1 (1.5) | 0.101 |
| cout | Mean (SD) | 583.2 (406.8) | 905.0 (945.5) | 0.102 |
| bonus | Mean (SD) | 55.8 (12.5) | 54.1 (10.1) | 0.550 |
| marque | Toyota | 7 (25.9) | 10 (26.3) | 0.054 |
| Mazda | 6 (22.2) | 17 (44.7) | ||
| Ford | 0 (0.0) | 2 (5.3) | ||
| Renault | 2 (7.4) | 1 (2.6) | ||
| Nissan | 0 (0.0) | 0 (0.0) | ||
| Peugeot | 0 (0.0) | 2 (5.3) | ||
| Honda | 0 (0.0) | 1 (2.6) | ||
| Bmw | 0 (0.0) | 1 (2.6) | ||
| Range | 10 (37.0) | 4 (10.5) | ||
| Volkswagen | 2 (7.4) | 0 (0.0) | ||
| Suziki | 0 (0.0) | 0 (0.0) | ||
| region | Abidjan | 2 (7.4) | 0 (0.0) | 0.148 |
| Bouake | 0 (0.0) | 0 (0.0) | ||
| Yakro | 0 (0.0) | 0 (0.0) | ||
| korhogo | 0 (0.0) | 1 (2.6) | ||
| Odienne | 0 (0.0) | 0 (0.0) | ||
| Man | 0 (0.0) | 1 (2.6) | ||
| Daloa | 1 (3.7) | 1 (2.6) | ||
| Duekoue | 0 (0.0) | 0 (0.0) | ||
| San pedro | 0 (0.0) | 0 (0.0) | ||
| Lahou | 0 (0.0) | 0 (0.0) | ||
| Soubre | 3 (11.1) | 0 (0.0) | ||
| Aboisso | 0 (0.0) | 1 (2.6) | ||
| Bongouanou | 0 (0.0) | 0 (0.0) | ||
| Daoukro | 21 (77.8) | 34 (89.5) | ||
| zone | Sud | 7 (25.9) | 13 (34.2) | 0.238 |
| Nord | 4 (14.8) | 5 (13.2) | ||
| Centre | 13 (48.1) | 11 (28.9) | ||
| Est | 3 (11.1) | 4 (10.5) | ||
| Cotière | 0 (0.0) | 5 (13.2) | ||
| Ouest | 0 (0.0) | 0 (0.0) |
)
res.m<-rstudent(modeLog)
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(modeLog)
sum(as.numeric(abs(res.m)<=3))/nrow(RegLogTrain)*100
## [1] 100