PROJET DE MODELISATION DE LA SINISTRALITE

I-Introduction

Dans ce projet, notre objectif est de modéliser le risque d’un portefeuille d’assurance. une fois que cette étape est franchie, il nous sera facile de faire la tarification, c’est-à-dire la détermination des primes.

Etant donné que le risque se décompose en deux dimensions: la fréquence de sinistres (le nombre de sinistres) et la sévérité (le coût moyen), il sera question de trouver les modèles qui optimisent ces derniers, c’est-dire des variables de la base de données qui les expliquent.

II-Importation et étude descriptive de la base de données

library(haven)

## Warning: le package 'haven' a été compilé avec la version R 4.1.1

base5 <- read_sas("~/Projet IFA/TP Assurance/base5.sas7bdat", 
    NULL)
head(base5,5)

Variables et leur signification

str(base5)

## tibble [16,082 x 19] (S3: tbl_df/tbl/data.frame)
##  $ NAP     : num [1:16082] 83 916 550 89 233 666 80 666 173 474 ...
##  $ PERMIS  : num [1:16082] 332 333 173 364 426 429 461 462 405 386 ...
##  $ DEB_IMAG: Date[1:16082], format: "2004-01-01" "2004-02-01" ...
##  $ FIN_IMAG: Date[1:16082], format: "2004-02-01" NA ...
##  $ SEX     : chr [1:16082] "M" "M" "M" "F" ...
##  $ STATUT  : chr [1:16082] "A" "A" "A" "A" ...
##  $ CSP     : num [1:16082] 50 50 50 55 60 60 48 48 50 55 ...
##  $ USAGE   : num [1:16082] 3 3 2 2 1 1 3 3 2 2 ...
##  $ AGECOND : num [1:16082] 46 46 32 52 57 57 58 58 53 57 ...
##  $ K8000   : chr [1:16082] "N" "N" "N" "N" ...
##  $ CHARGE  : num [1:16082] 0 0 0 0 0 0 0 0 0 0 ...
##  $ SINAP1  : num [1:16082] 0 0 0 0 0 0 0 0 0 0 ...
##  $ SINAP2  : num [1:16082] 1 1 2 0 0 0 1 1 0 0 ...
##  $ SINAP3  : num [1:16082] 0 0 0 0 0 0 0 0 0 0 ...
##  $ SINAP4  : num [1:16082] 0 0 0 0 0 0 0 0 0 0 ...
##  $ SINAP5  : num [1:16082] 0 0 0 0 0 0 2 2 0 1 ...
##  $ SINAP6  : num [1:16082] 0 0 0 0 0 0 1 1 0 0 ...
##  $ ZONE    : num [1:16082] 9 9 7 8 7 7 5 5 10 6 ...
##  $ RM      : num [1:16082] 50 50 68 50 50 50 50 50 50 50 ...

Notre base de données est un portefeuille d’assurance automobile composée de 19 variables dont 16082 observations.Une ligne correspond à un contrat.

age du conducteur: agecond, ancienneté de permis: permis, sexe du conducteur :sex, statut matrimonial :statut, catégorie socio-professionnelle :csp, usage du véhicule :usage, option kilométrage limité :k8000, zone géographique :zone, coefficient de réduction majoration (bonus/malus) :RM,date de début d’image :deb_imag,date de fin d’image :fin_imag, nombre d’années-police :nap, nombre de sinistres responsables dans les 4 années précédent l’image :sinap1, nombre de sinistres non responsables dans les 4 années précédent l’image :sinap2, nombre de sinistres parking dans les 4 années précédent l’image :sinap3, nombre de sinistres incendie/vol dans les 4 années précédent l’image :sinap4, nombre de sinistres bris de glace dans les 4 années précédent l’image sinap5, nombre de mises en demeure dans les 4 années précédent l’image :sinap6, charge de sinistres :charge.

Statistique de base

summary(base5)

##       NAP           PERMIS         DEB_IMAG             FIN_IMAG         
##  Min.   :   1   Min.   : 24.0   Min.   :2004-01-01   Min.   :2004-01-02  
##  1st Qu.: 173   1st Qu.:238.0   1st Qu.:2004-01-01   1st Qu.:2004-05-01  
##  Median : 408   Median :361.0   Median :2004-04-09   Median :2004-07-16  
##  Mean   : 425   Mean   :355.4   Mean   :2004-05-05   Mean   :2004-07-19  
##  3rd Qu.: 666   3rd Qu.:466.0   3rd Qu.:2004-08-05   3rd Qu.:2004-10-12  
##  Max.   :1000   Max.   :844.0   Max.   :2004-12-31   Max.   :2004-12-31  
##                                                      NA's   :7873        
##      SEX               STATUT               CSP            USAGE      
##  Length:16082       Length:16082       Min.   : 1.00   Min.   :1.000  
##  Class :character   Class :character   1st Qu.:50.00   1st Qu.:1.000  
##  Mode  :character   Mode  :character   Median :50.00   Median :2.000  
##                                        Mean   :49.92   Mean   :1.895  
##                                        3rd Qu.:60.00   3rd Qu.:2.000  
##                                        Max.   :66.00   Max.   :4.000  
##                                                                       
##     AGECOND         K8000               CHARGE            SINAP1      
##  Min.   :20.00   Length:16082       Min.   :-1842.0   Min.   :0.0000  
##  1st Qu.:40.00   Class :character   1st Qu.:    0.0   1st Qu.:0.0000  
##  Median :52.00   Mode  :character   Median :    0.0   Median :0.0000  
##  Mean   :50.95                      Mean   :  201.9   Mean   :0.2654  
##  3rd Qu.:61.00                      3rd Qu.:    0.0   3rd Qu.:0.0000  
##  Max.   :95.00                      Max.   :95151.0   Max.   :4.0000  
##                                                                       
##      SINAP2           SINAP3            SINAP4            SINAP5      
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.00000   Median :0.00000   Median :0.0000  
##  Mean   :0.3166   Mean   :0.08494   Mean   :0.07499   Mean   :0.4205  
##  3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:1.0000  
##  Max.   :7.0000   Max.   :4.00000   Max.   :3.00000   Max.   :8.0000  
##                                                                       
##      SINAP6            ZONE              RM        
##  Min.   :0.0000   Min.   : 1.000   Min.   : 50.00  
##  1st Qu.:0.0000   1st Qu.: 6.000   1st Qu.: 50.00  
##  Median :0.0000   Median : 7.000   Median : 50.00  
##  Mean   :0.2223   Mean   : 7.835   Mean   : 57.21  
##  3rd Qu.:0.0000   3rd Qu.:10.000   3rd Qu.: 58.00  
##  Max.   :5.0000   Max.   :13.000   Max.   :185.00  
## 

On remarque que la variable charge admet des valeurs négatives et nulles. Donc il est important de les enlever pour une bonne étude.

Connaissance des valeurs manquantes

library(funModeling)

## Warning: le package 'funModeling' a été compilé avec la version R 4.1.3

## Le chargement a nécessité le package : Hmisc

## Warning: le package 'Hmisc' a été compilé avec la version R 4.1.1

## Le chargement a nécessité le package : lattice

## Le chargement a nécessité le package : survival

## Le chargement a nécessité le package : Formula

## Warning: le package 'Formula' a été compilé avec la version R 4.1.1

## Le chargement a nécessité le package : ggplot2

## Warning: le package 'ggplot2' a été compilé avec la version R 4.1.2

## 
## Attachement du package : 'Hmisc'

## Les objets suivants sont masqués depuis 'package:base':
## 
##     format.pval, units

## funModeling v.1.9.4 :)
## Examples and tutorials at livebook.datascienceheroes.com
##  / Now in Spanish: librovivodecienciadedatos.ai

df_status(base5)

##    variable q_zeros p_zeros q_na  p_na q_inf p_inf      type unique
## 1       NAP       0    0.00    0  0.00     0     0   numeric    752
## 2    PERMIS       0    0.00    0  0.00     0     0   numeric    736
## 3  DEB_IMAG       0    0.00    0  0.00     0     0      Date    364
## 4  FIN_IMAG       0    0.00 7873 48.96     0     0      Date    364
## 5       SEX       0    0.00    0  0.00     0     0 character      2
## 6    STATUT       0    0.00    0  0.00     0     0 character      2
## 7       CSP       0    0.00    0  0.00     0     0   numeric     33
## 8     USAGE       0    0.00    0  0.00     0     0   numeric      4
## 9   AGECOND       0    0.00    0  0.00     0     0   numeric     73
## 10    K8000       0    0.00    0  0.00     0     0 character      2
## 11   CHARGE   14425   89.70    0  0.00     0     0   numeric   1425
## 12   SINAP1   12425   77.26    0  0.00     0     0   numeric      5
## 13   SINAP2   11969   74.42    0  0.00     0     0   numeric      8
## 14   SINAP3   14854   92.36    0  0.00     0     0   numeric      5
## 15   SINAP4   14973   93.10    0  0.00     0     0   numeric      4
## 16   SINAP5   10843   67.42    0  0.00     0     0   numeric      8
## 17   SINAP6   13773   85.64    0  0.00     0     0   numeric      6
## 18     ZONE       0    0.00    0  0.00     0     0   numeric     13
## 19       RM       0    0.00    0  0.00     0     0   numeric     75

library(naniar)
gg_miss_var(base5)

La variable FIN_IMG contient des valeurs manquantes, soit 48.96 %. On doit les enlever ou remplacer; mais on préfère les enlever pour ne peut pas corrompre les données.

Traitement de la base de données

base4=na.omit(base5)
attach(base4)
base4=subset(base4,CHARGE>0)
str(base4)

## tibble [760 x 19] (S3: tbl_df/tbl/data.frame)
##  $ NAP     : num [1:760] 761 23 833 957 482 998 97 707 250 250 ...
##  $ PERMIS  : num [1:760] 357 271 539 302 343 276 446 352 614 412 ...
##  $ DEB_IMAG: Date[1:760], format: "2004-01-01" "2004-10-23" ...
##  $ FIN_IMAG: Date[1:760], format: "2004-10-05" "2004-11-01" ...
##  $ SEX     : chr [1:760] "F" "M" "M" "F" ...
##  $ STATUT  : chr [1:760] "C" "A" "A" "A" ...
##  $ CSP     : num [1:760] 55 50 60 55 1 46 48 60 60 50 ...
##  $ USAGE   : num [1:760] 2 2 1 2 3 3 3 1 1 1 ...
##  $ AGECOND : num [1:760] 47 49 75 49 48 41 57 63 72 62 ...
##  $ K8000   : chr [1:760] "N" "N" "N" "N" ...
##  $ CHARGE  : num [1:760] 2018 357 645 1200 2667 ...
##  $ SINAP1  : num [1:760] 0 0 0 0 0 0 0 0 0 0 ...
##  $ SINAP2  : num [1:760] 1 1 0 0 0 0 1 0 0 0 ...
##  $ SINAP3  : num [1:760] 0 0 0 0 1 0 0 0 0 0 ...
##  $ SINAP4  : num [1:760] 0 0 0 0 0 0 0 0 0 0 ...
##  $ SINAP5  : num [1:760] 0 0 0 1 1 0 0 1 0 0 ...
##  $ SINAP6  : num [1:760] 0 2 0 0 0 0 0 1 0 0 ...
##  $ ZONE    : num [1:760] 6 8 9 8 9 6 11 4 9 10 ...
##  $ RM      : num [1:760] 50 50 50 54 64 60 50 54 50 52 ...
##  - attr(*, "na.action")= 'omit' Named int [1:7873] 2 4 6 8 9 11 14 15 18 20 ...
##   ..- attr(*, "names")= chr [1:7873] "2" "4" "6" "8" ...

Après traitement, notre base de données contient maintenant 760 observations, 19 variables.

Etude de la corrélation entre les variables quantitatives

library(corrplot)
Mat=base4[,c(1,2,7,8,18,19)]
corrplot(cor(Mat), type="upper", order="hclust", tl.col="black", tl.srt=45)

Il y a une coorélation disons moyenne mais négative entre CSP et USAGE.

library(FactoMineR)
library(factoextra)
Var_qual=base4[,c(5,6,10)]
Var_quali=FactoMineR::MCA(Var_qual,graph = FALSE)
factoextra::fviz_mca_var(Var_quali,repel = TRUE,choice = "mca.cor", ggtheme = theme_minimal ())

Plus les variables sont proches d’un axe, plus elles sont corrélées entre elles avec une dimension.

Etudes des variables

attach(base4)

sort(table(AGECOND))

## AGECOND
## 82 85 21 24 77 79 80 83 20 78 81 22 23 72 70 76 32 64 71 73 38 74 25 26 68 75 
##  1  1  2  2  2  2  2  2  3  3  3  5  5  5  6  6  7  7  7  7  8  8  9  9  9  9 
## 65 66 69 29 58 27 33 39 50 67 41 45 63 28 31 34 62 30 35 36 37 44 46 52 55 40 
## 10 10 10 11 11 13 13 13 13 13 14 14 14 15 15 15 15 16 16 16 16 16 16 16 16 17 
## 53 59 43 42 48 49 60 61 47 54 57 51 56 
## 17 17 18 19 19 19 19 21 22 22 22 24 27

Nombre=c( 1,  1 , 2  ,2 ,2  ,2 , 2 , 2  ,3 , 3 , 3 , 5 , 5,  5 , 6  ,6  ,7  ,7  ,7  ,7 , 8 , 8,  9 , 9  ,9  ,9 ,10 ,10 ,10 ,11,11, 13, 13 ,13, 13, 13 ,14 ,14, 14, 15 ,15 ,15, 15, 16 ,16 ,16, 16 ,16 ,16 ,16 ,16 ,17 ,17, 17 ,18, 19 ,19 ,19 ,19 ,21,22, 22, 22, 24 ,27 )

AGE=c( 82, 85, 21, 24, 77, 79 ,80, 83 ,20 ,78, 81 ,22 ,23 ,72, 70, 76 ,32, 64 ,71, 73 ,38 ,74, 25 ,26 ,68 ,75 ,65, 66 ,69, 29,58 ,27, 33 ,39, 50, 67, 41 ,45 ,63 ,28 ,31 ,34, 62, 30, 35, 36, 37 ,44, 46 ,52 ,55 ,40 ,53 ,59, 43 ,42 ,48 ,49 ,60 ,61 ,47, 54, 57 ,51, 56 )

library(plotly)

## Warning: le package 'plotly' a été compilé avec la version R 4.1.2

## 
## Attachement du package : 'plotly'

## L'objet suivant est masqué depuis 'package:Hmisc':
## 
##     subplot

## L'objet suivant est masqué depuis 'package:ggplot2':
## 
##     last_plot

## L'objet suivant est masqué depuis 'package:stats':
## 
##     filter

## L'objet suivant est masqué depuis 'package:graphics':
## 
##     layout

plot_ly(base5,x=~AGE,y=~Nombre,type="bar") %>% 
  layout(xaxis=list(title='AGE'),title="Répartition selon l'âge des conducteurs")

27 assurés ont 56 ans.Ce sont les majoritaires.

sort(table(STATUT))

## STATUT
##   C   A 
## 116 644

library(plotly)
effectif=c(644,116)
plot_ly(base5,x=~c("A",'C'),y=~effectif,type="bar") %>% 
  layout(title='Répartition selon le statut', xaxis=list(title='Statut'))

Le statut matrimonial dominant est A. 644 assurés sont de ce statut

sort(table(SEX))

## SEX
##   F   M 
## 264 496

Effectif=c(264 ,496)
sexe=c("Femmes",'Hommes')
plot_ly(base5,x=~sexe,y=~Effectif,type="bar") %>% 
  layout(title='Répartition selon le sexe')

La majorité des assurés sont des hommes, ils sont 496.

sort(table(ZONE))

## ZONE
##  12  13   2   3   8   4   5   6  11   9  10   7 
##   1   1   8  10  37  39  46  93  96 109 151 169

Effectif=c(1  , 1 ,  8 , 10 , 37 , 39 , 46  ,93 , 96 ,109 ,151 ,169 )
zone=c(12 , 13  , 2  , 3  , 8 ,  4 ,  5 ,  6 , 11 ,  9  ,10 ,  7 )
plot_ly(base5,x=~zone, y=~Effectif, type="bar") %>% 
  layout(title='Répartition selon la zone géographique',xaxis=list(title='zone géographique'))

La majorité des assurés sont dans la zone 7. Ils sont au nombre de 169.

sort(table(K8000))

## K8000
##   O   N 
##  43 717

Effectif=c(43,717)
plot_ly(base5,x=~c("o","N"), y=~Effectif, type="bar") %>% 
  layout(title='Répartition selon l\'option de kilométrage ', xaxis=list(title='option'))

L’option de kilométrage la plus répandue est N, 717 sur l’effectif.

sort(table(USAGE))

## USAGE
##   4   3   1   2 
##  43 152 217 348

Effectif=c(43, 152, 217 ,348 )
plot_ly(base5,x=~c(4  ,  3  ,  1,    2 ), y=~Effectif, type="bar") %>% 
  layout(title='Répartition selon l\'usage des véhicules ', xaxis=list(title='usage'))

L’usage 2 est la plus représentée.

sort(table(RM))

## RM
##  73  74  77  82  89  91  92  93 112 125 138 156  67  69 148  53  61  66  71  75 
##   1   1   1   1   1   1   1   1   1   1   1   1   2   2   2   3   3   3   3   3 
## 100 106  63  56  58  59 118  52  55  64  62  60  68  76  51  95  54  85  72  57 
##   3   3   4   5   5   5   5   6   6  11  12  13  13  13  15  15  17  17  20  21 
##  80  90  50 
##  24  34 460

effectif=c(1 ,  1 ,  1,   1  , 1  , 1   ,1 ,  1 ,  1 ,  1  , 1   ,1 ,  2  , 2  , 2  , 3  , 3 ,  3  , 3 ,  3  , 3  , 3,4 ,  5 ,  5,   5   ,5  , 6 ,  6 , 11 , 12,  13 , 13 , 13  ,15 , 15 , 17  ,17 , 20  ,21 , 24 , 34 ,460 )
rm=c( 73 , 74  ,77 , 82  ,89  ,91,  92 , 93 ,112 ,125 ,138 ,156 , 67 , 69 ,148 , 53 , 61 , 66,  71 , 75, 100 ,106,63  ,56 , 58  ,59, 118 , 52  ,55 , 64 , 62,  60  ,68  ,76  ,51  ,95,  54 , 85  ,72  ,57,  80 , 90 , 50 )

plot_ly(base5,x=~rm, y=~effectif, type="bar") %>% 
  layout(title='Répartition selon le coefficient de Réduction-Majoration', xaxis=list(title='RM'))

460 assurés sont bénéficiaires d’un coefficient de réduction-majoration égal à 50. C’est donc le RM dominant.

plot_ly(base4,type="histogram",x=~NAP) %>% 
  layout(title="Histogramme du nombre de police")

III-Modélisation de la fréquence de sinistre

NOMBRE=SINAP1+SINAP2+SINAP3+SINAP4+SINAP5+SINAP6
sort(table(NOMBRE))

## NOMBRE
##   9   7   6   5   4   3   2   0   1 
##   1   2   8  24  57  95 163 171 239

effectif=c(1 ,  2 ,  8  ,24,  57  ,95 ,163 ,171, 239  )
library(plotly)
plot_ly(type = 'bar',x=~c( 9   ,7 ,  6  , 5 ,  4  , 3  , 2 ,  0  , 1 ),y=~effectif) %>% 
  layout(title='Répartition du nombre de sinistres',
         xaxis=list(title='Nombre de sinistre'))

Nous avons en abscisse le nombre de sinistres et en ordonnée le nombre d’observations total associé au nombre de sinistres. Par exemple en 0, nous avons 171 assurés avec aucun sinistre. Il ya 239 qui ont un sinistre.

plot_ly(type = "pie",values=c(1 ,  2 ,  8  ,24,  57  ,95 ,163 ,171, 239  ),
        labels=c( 9   ,7 ,  6  , 5 ,  4  , 3  , 2 ,  0  , 1 ),texinfo="label+percent") 

## Warning: 'pie' objects don't have these attributes: 'texinfo'
## Valid attributes include:
## '_deprecated', 'automargin', 'customdata', 'customdatasrc', 'direction', 'dlabel', 'domain', 'hole', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'insidetextfont', 'insidetextorientation', 'label0', 'labels', 'labelssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'marker', 'meta', 'metasrc', 'name', 'opacity', 'outsidetextfont', 'pull', 'pullsrc', 'rotation', 'scalegroup', 'showlegend', 'sort', 'stream', 'text', 'textfont', 'textinfo', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'title', 'transforms', 'type', 'uid', 'uirevision', 'values', 'valuessrc', 'visible', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

22,5 % des assurés ont 0 sinistre.

Rep=xtabs(~SEX+NOMBRE,base4)
Rep

##    NOMBRE
## SEX   0   1   2   3   4   5   6   7   9
##   F  60  86  58  33  14  10   2   1   0
##   M 111 153 105  62  43  14   6   1   1

barplot(Rep,main = "Répartition du portefeuille selon le sexe en terme de sinistralité
",
xlab = "Nombre de sinistres",
col = c("red","green")
)
legend("topleft",
c("F","H"),
fill = c("red","green")
)

On constate que les hommes causent plus d’accidents que les femmes. Certainement parce qu’elles sont plus prudentes et roulent à faible vitesse; contrairement aux hommes.

Rep=xtabs(~STATUT+NOMBRE,base4)
Rep

##       NOMBRE
## STATUT   0   1   2   3   4   5   6   7   9
##      A 148 195 145  83  46  21   4   1   1
##      C  23  44  18  12  11   3   4   1   0

barplot(Rep,main = "Répartition du portefeuille selon le statut en terme de sinistralité
",
xlab = "Nombre de sinistres",
col = c("blue","green")
)
legend("topleft",
c("A","C"),
fill = c("blue","green")
)

p=xtabs(~K8000+NOMBRE,base4)
p

##      NOMBRE
## K8000   0   1   2   3   4   5   6   7   9
##     N 152 227 154  93  56  24   8   2   1
##     O  19  12   9   2   1   0   0   0   0

barplot(p,main = "Répartition du portefeuille selon l'option en terme de sinistralité
",
xlab = "Nombre de sinistres",
col = c("yellow","green")
)
legend("topleft",
c("N","O"),
fill = c("yellow","green")
)

Choix de la loi de la fréquence des sinistres: quelques test

x=rpois(length(NOMBRE),mean(NOMBRE))
ks.test(NOMBRE,x)

## Warning in ks.test(NOMBRE, x): p-value will be approximate in the presence of
## ties

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  NOMBRE and x
## D = 0.038158, p-value = 0.6376
## alternative hypothesis: two-sided

D’après le résultat du test de Kolmokorov(test d’ajustement), la fréquence de sinistres suit la même loi que celle de la variable x. Donc elle suit une distribution de poisson.

qqplot(NOMBRE,x,main="Test d'adéquation",pch=19,col='blue') 

On peut supposer que le nombre de sinistres suit une distribution de poisson

Test d’indépendance du nombre de sinistres avec les variables

IND=NOMBRE
IND[NOMBRE>0]=rep(1,length(IND[NOMBRE>0]))
chisq.test(AGECOND,IND)

## 
##  Pearson's Chi-squared test
## 
## data:  AGECOND and IND
## X-squared = 69.94, df = 64, p-value = 0.285

chisq.test(SEX,IND)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  SEX and IND
## X-squared = 0.00033285, df = 1, p-value = 0.9854

chisq.test(USAGE,IND)

## 
##  Pearson's Chi-squared test
## 
## data:  USAGE and IND
## X-squared = 12.011, df = 3, p-value = 0.007347

chisq.test(CSP,IND)

## 
##  Pearson's Chi-squared test
## 
## data:  CSP and IND
## X-squared = 23.316, df = 17, p-value = 0.1392

chisq.test(STATUT,IND)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  STATUT and IND
## X-squared = 0.3944, df = 1, p-value = 0.53

chisq.test(K8000,IND)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  K8000 and IND
## X-squared = 11.01, df = 1, p-value = 0.0009064

chisq.test(base4$RM,IND)

## 
##  Pearson's Chi-squared test
## 
## data:  base4$RM and IND
## X-squared = 39.057, df = 42, p-value = 0.6009

chisq.test(ZONE,IND)

## 
##  Pearson's Chi-squared test
## 
## data:  ZONE and IND
## X-squared = 6.6429, df = 11, p-value = 0.8272

chisq.test(NAP,IND)

## 
##  Pearson's Chi-squared test
## 
## data:  NAP and IND
## X-squared = 327.87, df = 343, p-value = 0.7125

chisq.test(PERMIS,IND)

## 
##  Pearson's Chi-squared test
## 
## data:  PERMIS and IND
## X-squared = 443.15, df = 420, p-value = 0.2097

Pour les variables AGECOND, SEX, CSP, ZONE,PERMIS, STATUT, RM et NAP, on accepte l’hypothèse Ho car leurs p-values sont supérieures à 0.05. Donc elles sont indépendantes du nombre de sinistres.

Courbes de fréquence de sinistres en fonction de quelques variables

n=max(base4$AGECOND)-min(base4$AGECOND)
 fage=seq(0,n) #fage :nb de sinistres correspondant à chaque age
 for(i in 1:n){
   x=19+i #on a mis 19 car i commence par 1
   y=NOMBRE[base4$AGECOND==x]
   fage[i]=sum(y)
 }
 w=seq(0,n)+20
 fage1=fage/sum(NOMBRE)
 fréquence=fage1
 library(plotly)
 plotly::plot_ly(x=~w,y=~fréquence,type='scatter', mode="lines",base4) %>% 
   layout(title='courbe de la fréquence de sinistre en fonction de lage', xaxis=list(title='âge'))

 n=max(base4$ZONE)-min(base4$ZONE)
 fzone=seq(0,n) #fage :nb de sinistres correspondant à chaque zone
 p=min(base4$ZONE)
 for(i in 1:n){
   x=p-1+i #on a mis 19 car i commence par 1
   y=NOMBRE[base4$ZONE==x]
   fzone[i]=sum(y)
 }
 w=seq(0,n)+p
 fzone1=fzone/sum(NOMBRE)
 #plot(w,fzone1,col="lightslateblue",type="l",lwd=2)
 library(plotly)
 fréquence=fzone1
 plotly::plot_ly(x=~w,y=~fréquence,type='scatter', mode="lines",base4) %>% 
   layout(title='courbe de la fréquence de sinistre en fonction de la zone', xaxis=list(title='zone'))

 n=max(base4$CSP)-min(base4$CSP)
 fCSP=seq(0,n) #fage :nb de sinistres correspondent à chaque catégorie
 p=min(base4$CSP)
 for(i in 1:n){
   x=p-1+i #on a mis 19 car i commence par 1
   y=NOMBRE[base4$CSP==x]
   fCSP[i]=sum(y)
 }
 w=seq(0,n)+p
 fCSP1=fCSP/sum(NOMBRE)

  library(plotly)
 fréquence=fCSP1
 plotly::plot_ly(x=~w,y=~fréquence,type='scatter', mode="lines",base4) %>% 
   layout(title='courbe de la fréquence de sinistre en fonction de la CSP', xaxis=list(title='CSP'))

Modèle linéaire généralisé

Frequence=glm(NOMBRE~PERMIS+AGECOND+ZONE+SEX+RM+USAGE+CSP+NAP+STATUT,family=poisson(link="log"),data = base4)
print(summary(Frequence))

## 
## Call:
## glm(formula = NOMBRE ~ PERMIS + AGECOND + ZONE + SEX + RM + USAGE + 
##     CSP + NAP + STATUT, family = poisson(link = "log"), data = base4)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2693  -0.8742  -0.2604   0.5851   3.6139  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.0867696  0.2959725   0.293  0.76939    
## PERMIS       0.0018742  0.0006208   3.019  0.00254 ** 
## AGECOND     -0.0186687  0.0068671  -2.719  0.00656 ** 
## ZONE        -0.0253005  0.0127151  -1.990  0.04661 *  
## SEXM        -0.0087647  0.0614144  -0.143  0.88652    
## RM           0.0087813  0.0017212   5.102 3.37e-07 ***
## USAGE        0.1814516  0.0394427   4.600 4.22e-06 ***
## CSP          0.0034267  0.0022504   1.523  0.12784    
## NAP         -0.0002852  0.0001142  -2.497  0.01254 *  
## STATUTC     -0.0184665  0.0861354  -0.214  0.83024    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1061.45  on 759  degrees of freedom
## Residual deviance:  990.66  on 750  degrees of freedom
## AIC: 2509.5
## 
## Number of Fisher Scoring iterations: 5

Sélection des variables par la méthode backward

library(MASS)

## 
## Attachement du package : 'MASS'

## L'objet suivant est masqué depuis 'package:plotly':
## 
##     select

Frequence_sinistre=stepAIC(Frequence,direction = 'backward',k=log(nrow(base4)))

## Start:  AIC=2555.86
## NOMBRE ~ PERMIS + AGECOND + ZONE + SEX + RM + USAGE + CSP + NAP + 
##     STATUT
## 
##           Df Deviance    AIC
## - SEX      1   990.68 2549.2
## - STATUT   1   990.71 2549.3
## - CSP      1   993.03 2551.6
## - ZONE     1   994.60 2553.2
## - NAP      1   996.88 2555.4
## <none>         990.66 2555.9
## - AGECOND  1   998.54 2557.1
## - PERMIS   1  1000.33 2558.9
## - USAGE    1  1011.29 2569.9
## - RM       1  1013.93 2572.5
## 
## Step:  AIC=2549.25
## NOMBRE ~ PERMIS + AGECOND + ZONE + RM + USAGE + CSP + NAP + STATUT
## 
##           Df Deviance    AIC
## - STATUT   1   990.72 2542.7
## - CSP      1   993.03 2545.0
## - ZONE     1   994.62 2546.6
## - NAP      1   996.89 2548.8
## <none>         990.68 2549.2
## - AGECOND  1   998.61 2550.6
## - PERMIS   1  1000.56 2552.5
## - USAGE    1  1011.42 2563.4
## - RM       1  1013.93 2565.9
## 
## Step:  AIC=2542.66
## NOMBRE ~ PERMIS + AGECOND + ZONE + RM + USAGE + CSP + NAP
## 
##           Df Deviance    AIC
## - CSP      1   993.08 2538.4
## - ZONE     1   994.68 2540.0
## - NAP      1   996.92 2542.2
## <none>         990.72 2542.7
## - AGECOND  1   998.67 2544.0
## - PERMIS   1  1000.58 2545.9
## - USAGE    1  1012.12 2557.4
## - RM       1  1014.30 2559.6
## 
## Step:  AIC=2538.38
## NOMBRE ~ PERMIS + AGECOND + ZONE + RM + USAGE + NAP
## 
##           Df Deviance    AIC
## - ZONE     1   997.67 2536.3
## - NAP      1   999.20 2537.9
## <none>         993.08 2538.4
## - AGECOND  1  1000.51 2539.2
## - PERMIS   1  1002.56 2541.2
## - USAGE    1  1012.86 2551.5
## - RM       1  1017.03 2555.7
## 
## Step:  AIC=2536.34
## NOMBRE ~ PERMIS + AGECOND + RM + USAGE + NAP
## 
##           Df Deviance    AIC
## - NAP      1  1003.67 2535.7
## <none>         997.67 2536.3
## - AGECOND  1  1005.11 2537.2
## - PERMIS   1  1007.18 2539.2
## - USAGE    1  1018.65 2550.7
## - RM       1  1020.38 2552.4
## 
## Step:  AIC=2535.7
## NOMBRE ~ PERMIS + AGECOND + RM + USAGE
## 
##           Df Deviance    AIC
## <none>         1003.7 2535.7
## - AGECOND  1   1011.3 2536.7
## - PERMIS   1   1013.3 2538.7
## - USAGE    1   1024.6 2550.0
## - RM       1   1027.6 2553.0

print(summary(Frequence_sinistre))

## 
## Call:
## glm(formula = NOMBRE ~ PERMIS + AGECOND + RM + USAGE, family = poisson(link = "log"), 
##     data = base4)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3144  -0.8617  -0.2934   0.5698   3.8554  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.0841171  0.2284695  -0.368  0.71274    
## PERMIS       0.0018387  0.0006107   3.011  0.00260 ** 
## AGECOND     -0.0179794  0.0067151  -2.677  0.00742 ** 
## RM           0.0088164  0.0017029   5.177 2.25e-07 ***
## USAGE        0.1571377  0.0339235   4.632 3.62e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1061.5  on 759  degrees of freedom
## Residual deviance: 1003.7  on 755  degrees of freedom
## AIC: 2512.5
## 
## Number of Fisher Scoring iterations: 5

D’âprès ce qui précède, les variables qui optimisent notre modèle sont: PERMIS, AGECOND, RM et USAGE.

Frequence_sinistre = a0 + a1.PERMIS + a2.AGECOND + a3.RM + a4.USAGE

Validation du modèle

La validation d’un modèle GLM se fait par l’analyse des résidus du modèle GLM. Les résidus permettant de connaitre la qualité d’ajustement d’un modèle sont les résidus de déviance et de Pearson.

valeur=residuals.glm(Frequence_sinistre,type = 'deviance')
plot(valeur)

valeur=residuals.glm(Frequence_sinistre,type = 'pearson')
plot(valeur)

Une bonne qualité d’ajustement est traduite par l’appartenance des résidus de deviance standardisés à l’intervalle [−2 : 2]. Sur ces deux graphiques, nous observons que la plupart des points appartiennent à cet intervalle. Ce qui confirme que notre modèle est juste.

Prédiction

predict(Frequence_sinistre,newdata=data.frame(PERMIS=600,AGECOND=25,RM=80,USAGE=2),
        se.fit = TRUE,type='response',level=0.95,interval='prediction')

## $fit
##        1 
## 4.900012 
## 
## $se.fit
##        1 
## 1.583007 
## 
## $residual.scale
## [1] 1

Pour PERMIS=600,AGECOND=25,RM=80,USAGE=2, on aurait à peu près 5 sinistres.

Classification

library(rpart)
library(rpart.plot)
tre=rpart(NOMBRE~PERMIS+AGECOND+ZONE+SEX+base4$RM+USAGE+CSP+NAP+STATUT)
tre

## n= 760 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
## 1) root 760 1606.10400 1.669737  
##   2) base4$RM< 54.5 501  941.24950 1.500998  
##     4) USAGE< 2.5 367  561.42780 1.302452 *
##     5) USAGE>=2.5 134  325.73130 2.044776 *
##   3) base4$RM>=54.5 259  622.99610 1.996139  
##     6) PERMIS< 414.5 220  512.92730 1.881818 *
##     7) PERMIS>=414.5 39   90.97436 2.641026 *

rpart.plot(tre)

IV-Modélisation du coût de sinistre

Voici quelques lois usuelles pour la sévérité

z<- seq(0,5,.01)
y1=dlnorm(z,-0.5,1)
y2=dgamma(z,2,2)
y3=dweibull(z,2,2/sqrt(pi))
y4=dexp(z,1)
y5=dnorm(z,0,1)

library(plotly)
img=plot_ly(x=~z,y=~y1,type = "scatter",mode="lines",name = "loi log-normale")
  img=img %>%  add_trace(y=~y2,name = "loi de gamma")
  img=img %>%  add_trace(y=~y3,name = "loi de weibull")
  img=img %>%  add_trace(y=~y4,name = "loi exponentielle")
  img=img %>% add_trace(y=~y5, name='Loi normale')  
  img=img %>% layout(title="Densité de lois continues")  
img

On suppose la charge suit la loi de Gamma.

Ma_charge=glm(CHARGE~PERMIS+AGECOND+ZONE+SEX+RM+USAGE+CSP+NAP+STATUT,family=Gamma(link='log'),base4)

print(summary(Ma_charge))

## 
## Call:
## glm(formula = CHARGE ~ PERMIS + AGECOND + ZONE + SEX + RM + USAGE + 
##     CSP + NAP + STATUT, family = Gamma(link = "log"), data = base4)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.1148  -1.4780  -0.8755  -0.1299   5.6108  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.2005866  0.8129769   7.627 7.29e-14 ***
## PERMIS       0.0009441  0.0014631   0.645 0.518958    
## AGECOND      0.0077504  0.0158360   0.489 0.624691    
## ZONE         0.0263194  0.0336215   0.783 0.433985    
## SEXM         0.0172541  0.1602230   0.108 0.914272    
## RM           0.0193209  0.0053763   3.594 0.000347 ***
## USAGE        0.0464275  0.1088474   0.427 0.669839    
## CSP         -0.0064873  0.0061476  -1.055 0.291643    
## NAP         -0.0004816  0.0003018  -1.596 0.110996    
## STATUTC      0.3176955  0.2301449   1.380 0.167870    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 4.148662)
## 
##     Null deviance: 1665.5  on 759  degrees of freedom
## Residual deviance: 1562.8  on 750  degrees of freedom
## AIC: 13301
## 
## Number of Fisher Scoring iterations: 10

library(MASS)
cout_moyen=stepAIC(Ma_charge,direction = 'backward',k=log(nrow(base4)))

## Start:  AIC=13347.75
## CHARGE ~ PERMIS + AGECOND + ZONE + SEX + RM + USAGE + CSP + NAP + 
##     STATUT
## 
##           Df Deviance   AIC
## - SEX      1   1562.8 13341
## - USAGE    1   1563.6 13341
## - AGECOND  1   1564.0 13341
## - PERMIS   1   1564.6 13342
## - ZONE     1   1565.0 13342
## - CSP      1   1567.6 13342
## - STATUT   1   1570.2 13343
## - NAP      1   1571.9 13343
## <none>         1562.8 13348
## - RM       1   1623.0 13356
## 
## Step:  AIC=13341.15
## CHARGE ~ PERMIS + AGECOND + ZONE + RM + USAGE + CSP + NAP + STATUT
## 
##           Df Deviance   AIC
## - USAGE    1   1563.7 13335
## - AGECOND  1   1564.0 13335
## - PERMIS   1   1564.9 13335
## - ZONE     1   1565.1 13335
## - CSP      1   1567.6 13336
## - STATUT   1   1570.4 13336
## - NAP      1   1572.0 13337
## <none>         1562.8 13341
## - RM       1   1623.2 13349
## 
## Step:  AIC=13335.07
## CHARGE ~ PERMIS + AGECOND + ZONE + RM + CSP + NAP + STATUT
## 
##           Df Deviance   AIC
## - AGECOND  1   1564.6 13329
## - ZONE     1   1565.7 13329
## - PERMIS   1   1566.1 13329
## - STATUT   1   1570.7 13330
## - CSP      1   1573.0 13331
## - NAP      1   1573.0 13331
## <none>         1563.7 13335
## - RM       1   1624.8 13343
## 
## Step:  AIC=13328.98
## CHARGE ~ PERMIS + ZONE + RM + CSP + NAP + STATUT
## 
##          Df Deviance   AIC
## - ZONE    1   1567.0 13323
## - STATUT  1   1570.9 13324
## - CSP     1   1573.0 13324
## - NAP     1   1573.6 13324
## <none>        1564.6 13329
## - PERMIS  1   1594.9 13330
## - RM      1   1633.4 13339
## 
## Step:  AIC=13323.81
## CHARGE ~ PERMIS + RM + CSP + NAP + STATUT
## 
##          Df Deviance   AIC
## - STATUT  1   1573.5 13319
## - NAP     1   1576.0 13319
## - CSP     1   1576.1 13319
## <none>        1567.0 13324
## - PERMIS  1   1596.8 13324
## - RM      1   1641.7 13336
## 
## Step:  AIC=13321.21
## CHARGE ~ PERMIS + RM + CSP + NAP
## 
##          Df Deviance   AIC
## - CSP     1   1580.7 13316
## - NAP     1   1585.5 13318
## - PERMIS  1   1597.0 13320
## <none>        1573.5 13321
## - RM      1   1650.4 13334
## 
## Step:  AIC=13319.03
## CHARGE ~ PERMIS + RM + NAP
## 
##          Df Deviance   AIC
## - NAP     1   1593.4 13315
## - PERMIS  1   1600.3 13317
## <none>        1580.7 13319
## - RM      1   1652.4 13329
## 
## Step:  AIC=13320.19
## CHARGE ~ PERMIS + RM

print(summary(cout_moyen))

## 
## Call:
## glm(formula = CHARGE ~ PERMIS + RM, family = Gamma(link = "log"), 
##     data = base4)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.1009  -1.4732  -0.8808  -0.1716   6.6289  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.3070764  0.4453300  14.163  < 2e-16 ***
## PERMIS      0.0010479  0.0005584   1.877 0.060944 .  
## RM          0.0200815  0.0053481   3.755 0.000187 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 4.332298)
## 
##     Null deviance: 1665.5  on 759  degrees of freedom
## Residual deviance: 1593.4  on 757  degrees of freedom
## AIC: 13306
## 
## Number of Fisher Scoring iterations: 7

Les variables qui optimisent le coût de sinistre sont: PERMIS et RM cout_moyen =b0+ b1.PERMIS + b2.RM

Validation du modèle

valeur=residuals.glm(cout_moyen,type = 'deviance')
plot(valeur)

valeur=residuals.glm(cout_moyen,type = 'pearson')
plot(valeur)

Une bonne qualité d’ajustement est traduite par l’appartenance des résidus de deviance standardisés à l’intervalle [−2 : 2]. Sur ces deux graphiques, nous observons que la plupart des points appartiennent à cet intervalle. Ce qui confirme que notre modèle est juste.

Prédiction

predict(cout_moyen,newdata=data.frame(PERMIS=600,RM=80),
        se.fit = TRUE,type='response',level=0.95,interval='prediction')

## $fit
##        1 
## 5127.188 
## 
## $se.fit
##        1 
## 1207.907 
## 
## $residual.scale
## [1] 2.081417

Pour PERMIS=600,RM=80, on aurait une charge de 5127,188

Classification

library(rpart)
library(rpart.plot)
tre=rpart( CHARGE~PERMIS+AGECOND+ZONE+SEX+base4$RM+USAGE+CSP+NAP+STATUT)
tre

## n= 760 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 760 28901890000  2685.315  
##    2) base4$RM< 115 750 19665340000  2542.252  
##      4) AGECOND< 79.5 742 17860370000  2472.275  
##        8) NAP>=638 296  3253137000  1825.457 *
##        9) NAP< 638 446 14401210000  2901.552  
##         18) base4$RM< 55.5 289  7052572000  2294.390  
##           36) CSP>=13 269  3347903000  2007.688 *
##           37) CSP< 13 20  3385161000  6150.533  
##             74) NAP< 407.5 10    22438970  1625.171 *
##             75) NAP>=407.5 10  2953144000 10675.900 *
##         19) base4$RM>=55.5 157  7045983000  4019.196  
##           38) AGECOND< 63.5 145  4841357000  3569.214 *
##           39) AGECOND>=63.5 12  1820497000  9456.479 *
##      5) AGECOND>=79.5 8  1464328000  9032.686 *
##    3) base4$RM>=115 10  8069927000 13415.050 *

rpart.plot(tre)

V-Conclusion

Connaissant les modèles qui optimisent la fréquence et le coût de sinistres, il revient facile de calculer les primes en utilisant la prédiction des nouvelles valeurs, des formules du cours et les chargements nécessaires.