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

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

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

Traitement 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(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)

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(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
}

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(RegLog)

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(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)
}

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

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.

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.

Statistique sommaire la variable cible

freq(RegLog$accident)
##        n    % val%
## non 1314 47.5 47.5
## Oui 1451 52.5 52.5

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(RegLog.kNN)
# Echantillon d'apprentissage
RegLogTrain<- df[-c(2701:2765), ]
# Echantillon de test, prediction
RegLogTest<- 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)
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

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.

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

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.

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

Repré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(modeLog)

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(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"))

Matrice de confusion

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

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

Analyse des résidus

)

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

CONCLUSION