Introduction

Les données

Nous disposons de deux jeux de données initialement. freMTPLfreq de longueur 413169 qui répertorie :

summary(Contrats)
##     PolicyID         ClaimNb           Exposure            Power      
##  1      :     1   Min.   :0.00000   Min.   :0.002732   f      :95718  
##  2      :     1   1st Qu.:0.00000   1st Qu.:0.200000   g      :91198  
##  3      :     1   Median :0.00000   Median :0.540000   e      :77022  
##  4      :     1   Mean   :0.03916   Mean   :0.561088   d      :68014  
##  5      :     1   3rd Qu.:0.00000   3rd Qu.:1.000000   h      :26698  
##  6      :     1   Max.   :4.00000   Max.   :1.990000   j      :18038  
##  (Other):413163                                        (Other):36481  
##      CarAge          DriverAge    
##  Min.   :  0.000   Min.   :18.00  
##  1st Qu.:  3.000   1st Qu.:34.00  
##  Median :  7.000   Median :44.00  
##  Mean   :  7.532   Mean   :45.32  
##  3rd Qu.: 12.000   3rd Qu.:54.00  
##  Max.   :100.000   Max.   :99.00  
##                                   
##                                 Brand             Gas        
##  Fiat                              : 16723   Diesel :205945  
##  Japanese (except Nissan) or Korean: 79060   Regular:207224  
##  Mercedes, Chrysler or BMW         : 19280                   
##  Opel, General Motors or Ford      : 37402                   
##  other                             :  9866                   
##  Renault, Nissan or Citroen        :218200                   
##  Volkswagen, Audi, Skoda or Seat   : 32638                   
##      Region          Density     
##  R24    :160601   Min.   :    2  
##  R11    : 69791   1st Qu.:   67  
##  R53    : 42122   Median :  287  
##  R52    : 38751   Mean   : 1985  
##  R72    : 31329   3rd Qu.: 1410  
##  R31    : 27285   Max.   :27000  
##  (Other): 43290

L’autre jeu de données freMTPLsev de longueur 16181récapitule les montants de sinistres ClaimAmount pour chaque PolicyID.

summary(Sinistres)
##     PolicyID      ClaimAmount     
##  18652  :    4   Min.   :      2  
##  226856 :    4   1st Qu.:    698  
##  311457 :    4   Median :   1156  
##  15608  :    3   Mean   :   2130  
##  25363  :    3   3rd Qu.:   1243  
##  25580  :    3   Max.   :2036833  
##  (Other):16160

L’on constate que plusieurs lignes correspondent au même PolicyID. En effet, dans ce jeu de données, les montants des sinistres sont répertoriés un par un et non aggrégés par numéro de police. Nous allons donc regrouper les montants de sinistres par numéro de police. Nous utiliserons les librairies dplyr et tidyr pour une gestion des données plus efficace.

Sinistres = summarise(group_by(Sinistres,PolicyID),ClaimAmount =sum(ClaimAmount))

summary(Sinistres)
##     PolicyID      ClaimAmount     
##  33     :    1   Min.   :      2  
##  41     :    1   1st Qu.:    769  
##  92     :    1   Median :   1163  
##  96     :    1   Mean   :   2239  
##  142    :    1   3rd Qu.:   1279  
##  182    :    1   Max.   :2036833  
##  (Other):15384

Aprés avoir éliminé les doublons, nous avons maintenant 15390 lignes qui correspond au nombres d’assurés sinistrés.

Nous allons maintenant fusionner Sinistres et COntrats puis reconstituer deux nouveaux de jeux données, mydata.amount et mydata.freq ayant comme réponses respectives les montants de sinistres et les nombres de sinitres. Ces tables contiennent aussi toutes les autres variables explicatives à notre disposition.

### Fusion des deux jeux de données
mydata = merge(Contrats,Sinistres,by.x = "PolicyID", by.y = "PolicyID", all.x = TRUE)
mydata = replace(mydata, is.na(mydata), 0)
###Séparation en une base pour les montants de sinistres et une autre pour le nombre
mydata.amount =mydata %>%
  filter(ClaimAmount != 0) %>%
  select(- c(ClaimNb,PolicyID))
#summary(mydata.amount)
mydata.freq = mydata  %>%
  select(-c(ClaimAmount))
#summary(mydata.freq)

Nous procédons par la suite à la séparation des sinistres en graves et attritionnels. En effectuant un Mean Excess Plot, nous avons décidé de fixer un seuil à 20000. De ce fait il y a 112 sinistres qui seront considèrés comme étant graves et 99,31 % des sinistres sont considérés comme atritionnels.

### Choix d'un seuil pour les sinistres graves et séparation des données
meplot(mydata$ClaimAmount, col = "steelblue")
abline(v = 20000,col = "red")

#findthresh(mydata.amount$ClaimAmount, n = 133)# 133 sinistres considérés graves
#print(quantile(mydata.amount$ClaimAmount, probs = 0.99161))# proba sinistres attritionnels
seuil = 20000

Une fois notre seuil trouvé, nous séparons donc nos données. Nous précisons que cette séparation concerne que notre jeu de données mydata.amount qui nous servira à modéliser le coût des sinistres mais aussi mydata.freq.

policy.grave = mydata[mydata$ClaimAmount >= seuil,1]

mydata.amount.attritionnel = mydata.amount %>%
  filter(ClaimAmount <= seuil)
mydata.amount.grave = mydata.amount %>%
  filter(ClaimAmount > seuil)
mydata.freq = mydata.freq %>%
  filter(!(PolicyID %in% policy.grave)) %>%
  select(-PolicyID)

Avant de commencer la modélisation du coût et de la fréquence, nous allons d’abord verifier les corrélations entre les variables quantitatives pour éviter toute redondance entre celles-ci.

### Corrélation entre nos variables explicatives quantitatives

mydata.numeric = mydata %>%
  select(c(Exposure,CarAge,DriverAge,Density))

mat.cor = cor(mydata.numeric)

col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(mat.cor, method="color", col=col(200),  type="upper", addCoef.col = "black", 
         tl.col="black", tl.srt=45,diag=FALSE )

Les corrélations entre les variables explicactives quantitatives sont faibles, on les garde toutes donc. Nous allons commencer maintenant le calcul de la prime pure par le modèle fréquence-coût.

I. Modélisation par les arbres de régression et de classification

DAns cette partie, nous détaillerons entièrement le procédé qu’on utilise pour modéliser la fréquence et le montant des sinistres. Les autres parties seront donc plus succinctes. ### I.1 Modélisation de la fréquence Le jeu de données correspondant est mydata.freq. Nous le décomposerons en deux pour avoir un échantillon d’apprentissage et un autre de test.

set.seed(1234)
N = nrow(mydata.freq)

e = sort(sample(N, floor(N/2)))
e.1 = sort(sample(N, floor(N/5)))
freq.train = mydata.freq[e,]
freq.test = mydata.freq[-e,]
freq.train.1 = mydata.freq[e.1,]
freq.test.1 = mydata.freq[-e.1,]
summary(freq.train.1)
##     ClaimNb           Exposure            Power           CarAge       
##  Min.   :0.00000   Min.   :0.002732   f      :19174   Min.   :  0.000  
##  1st Qu.:0.00000   1st Qu.:0.200000   g      :18418   1st Qu.:  3.000  
##  Median :0.00000   Median :0.540000   e      :15235   Median :  7.000  
##  Mean   :0.04013   Mean   :0.563072   d      :13447   Mean   :  7.525  
##  3rd Qu.:0.00000   3rd Qu.:1.000000   h      : 5263   3rd Qu.: 12.000  
##  Max.   :4.00000   Max.   :1.750000   j      : 3658   Max.   :100.000  
##                                       (Other): 7412                    
##    DriverAge                                    Brand      
##  Min.   :18.00   Fiat                              : 3278  
##  1st Qu.:34.00   Japanese (except Nissan) or Korean:15779  
##  Median :44.00   Mercedes, Chrysler or BMW         : 3908  
##  Mean   :45.37   Opel, General Motors or Ford      : 7494  
##  3rd Qu.:54.00   other                             : 1982  
##  Max.   :99.00   Renault, Nissan or Citroen        :43718  
##                  Volkswagen, Audi, Skoda or Seat   : 6448  
##       Gas            Region         Density     
##  Diesel :41115   R24    :32077   Min.   :    2  
##  Regular:41492   R11    :13766   1st Qu.:   67  
##                  R53    : 8440   Median :  287  
##                  R52    : 7811   Mean   : 1969  
##                  R72    : 6285   3rd Qu.: 1414  
##                  R31    : 5515   Max.   :27000  
##                  (Other): 8713

Nous vérifions ensuite que nos échantillons reprennent bien les caractéristiques essentielles des données initiales.

A = matrix(c(mean(mydata.freq$ClaimNb), mean(freq.train$ClaimNb),mean(freq.test$ClaimNb),
             min(mydata.freq$ClaimNb),min(freq.train$ClaimNb),min(freq.test$ClaimNb),
             max(mydata.freq$ClaimNb),max(freq.train$ClaimNb),max(freq.test$ClaimNb)), nrow = 3)
rownames(A) = c("données initiales", "apprentissage", "test")
colnames(A) = c("moyenne", "min", "max")
A
##                      moyenne min max
## données initiales 0.03879816   0   4
## apprentissage     0.03881036   0   4
## test              0.03878597   0   4
A.1 = matrix(c(mean(mydata.freq$ClaimNb), mean(freq.train.1$ClaimNb),mean(freq.test.1$ClaimNb),
             min(mydata.freq$ClaimNb),min(freq.train.1$ClaimNb),min(freq.test.1$ClaimNb),
             max(mydata.freq$ClaimNb),max(freq.train.1$ClaimNb),max(freq.test.1$ClaimNb)), nrow = 3)
rownames(A.1) = c("données initiales", "apprentissage", "test")
colnames(A.1) = c("moyenne", "min", "max")
A.1
##                      moyenne min max
## données initiales 0.03879816   0   4
## apprentissage     0.04012977   0   4
## test              0.03846526   0   4

On peut procéder maintenant à l’implémentation de l’algorithme CART grâce à rpart().

set.seed(12345)

freq.cart = rpart(ClaimNb ~ . ,data =  freq.train, minbucket = 1000, cp = 0, weights = Exposure , method = "anova")
#summary(freq.cart)
fancyRpartPlot(freq.cart)
## Warning: labs do not fit even at cex 0.15, there may be some overplotting

#plot(nb.cart, uniform = TRUE)

Nous procédons à l’élagage de l’arbre pour obtenir l’arbre optimal.

freq.cart.cp.min = freq.cart.cp[freq.cart.cp[,4] == min(freq.cart.cp[,4]),1]
plot(freq.cart.cp[,1], freq.cart.cp[,4], type = "l")

La complexité cp optimal 1.563520910^{-4} est . L’élagage de l’arbre donne :

freq.cart.pruned = prune(freq.cart, cp = freq.cart.cp.min)
fancyRpartPlot(freq.cart.pruned)

#plot(nb.cart.pruned, uniform = TRUE)

Nous pouvons maintenant tester notre arbre sur l’échantillon de test.

freq.cart.test = predict(freq.cart.pruned, newdata = freq.test, type ="vector" )

Nous mesurons la qualité de la prédiction par la least square error :

lse = mean((freq.cart.test-freq.test$ClaimNb)^2)
lse
## [1] 0.04097941
mean(freq.cart.test)
## [1] 0.04055999

La modélisation avec les autres algorithmes nous permettra de mieux apprécier les mse obtenues dans ce chapitre.

II.2 Modélisation du montant des sinistres

Nous essayons à présent de modéliser le montant des sinistres.Nous séparons notre jeu de données mydata.amount.attritionnel en échantillons d’apprentissage amount.train (2/3 des données) et de test amount.test.

### Séparation apprentissage et test
L = nrow(mydata.amount.attritionnel)
set.seed(12345)
d = sort(sample(L, L/2))

amount.train = mydata.amount.attritionnel[d,]
amount.test = mydata.amount.attritionnel[-d,]

Nous vérifions ensuite que nos échantillons reprennent bien les caractéristiques essentielles des données initiales.

#Test sur les  échantillons
B = matrix(c(mean(mydata.amount.attritionnel$ClaimAmount), mean(amount.train$ClaimAmount),mean(amount.test$ClaimAmount),
             min(mydata.amount.attritionnel$ClaimAmount),min(amount.train$ClaimAmount),min(amount.test$ClaimAmount),
             max(mydata.amount.attritionnel$ClaimAmount),max(amount.train$ClaimAmount),max(amount.test$ClaimAmount)), nrow = 3)
rownames(B) = c("données initiales", "apprentissage", "test")
colnames(B) = c("moyenne", "min", "max")
B
##                    moyenne min   max
## données initiales 1455.986   2 19912
## apprentissage     1458.601   2 19912
## test              1453.370   2 19618

Arbre de régression pour le montant des sinistres.

set.seed(12345)

amount.cart = rpart(ClaimAmount ~. , data = amount.train ,  minbucket = 200, cp = 0, weights = Exposure )
fancyRpartPlot(amount.cart) 
## Warning: labs do not fit even at cex 0.15, there may be some overplotting

Comme pour la fréquence nous allons essayer d’optimiser l’arbre:

amount.cart.cp = printcp(amount.cart)
## 
## Regression tree:
## rpart(formula = ClaimAmount ~ ., data = amount.train, weights = Exposure, 
##     minbucket = 200, cp = 0)
## 
## Variables actually used in tree construction:
## [1] Brand     CarAge    Density   DriverAge Exposure  Power     Region   
## 
## Root node error: 1.353e+10/7628 = 1773668
## 
## n= 7628 
## 
##            CP nsplit rel error  xerror     xstd
## 1  0.00706769      0   1.00000 1.00098 0.086154
## 2  0.00369175      1   0.99293 0.99529 0.085600
## 3  0.00191862      2   0.98924 0.99345 0.085554
## 4  0.00176082      3   0.98732 0.99538 0.085433
## 5  0.00145701      4   0.98556 0.99648 0.085491
## 6  0.00137702      7   0.98119 0.99703 0.085456
## 7  0.00097413      8   0.97981 0.99950 0.085580
## 8  0.00094985      9   0.97884 1.00039 0.085414
## 9  0.00094805     10   0.97789 1.00042 0.085397
## 10 0.00092832     12   0.97599 1.00042 0.085397
## 11 0.00087747     13   0.97506 1.00096 0.085351
## 12 0.00084374     14   0.97419 1.00127 0.085285
## 13 0.00053099     15   0.97334 1.00123 0.085332
## 14 0.00042017     16   0.97281 1.00394 0.085479
## 15 0.00030606     17   0.97239 1.00399 0.085470
## 16 0.00027176     18   0.97209 1.00465 0.085579
## 17 0.00021782     19   0.97181 1.00471 0.085597
## 18 0.00000000     20   0.97160 1.00449 0.085563
amount.cart.cp.min = amount.cart.cp[amount.cart.cp[,4] == min(amount.cart.cp[,4]),1]
plot(amount.cart.cp[,1], amount.cart.cp[,4], type = "l")

En essayant d’élaguer l’abre, nous obtenons trop peu de feuilles. Nous allons donc considérer l’arbre optimal. De plus, nous verrons que la mse est meilleure avec l’abre maximal.

amount.cart.pruned = prune(amount.cart, cp = amount.cart.cp.min)

Prédiction et analyse des résidus sur l’échantillon de test.

amount.cart.pred = predict(amount.cart,amount.test)
amount.cart.pred.pruned = predict(amount.cart.pruned, amount.test)
summary(amount.cart.pred)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   991.6  1303.0  1360.0  1458.0  1677.0  2017.0
summary(amount.cart.pred.pruned)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1119    1367    1562    1428    1562    1562

Nous mesurons la qualité de la prédiction par la least square error :

lse.amount.cart = mean((amount.cart.pred-amount.test$ClaimAmount)^2)
lse.amount.cart.pruned = mean((amount.cart.pred.pruned-amount.test$ClaimAmount)^2)
lse.amount.cart
## [1] 3124853
lse.amount.cart.pruned
## [1] 3111924

Prime pure sans graves

sum(predict(freq.cart,mydata)*predict(amount.cart,mydata))
## [1] 23658419

III. Modélisation par les forêts aléatoires (rf)

Nous utilserons dans cette partie les échantillons d’apprentissage et de test créés dans le chapitre sur les arbres de régression et de classification.

Pour les fores aléatoires, nous utiliserons le package éponyme randomForest.

III.1 Modélisation de la fréquence

III.2 Modélisation du montant des sinistres

Apprentissage

set.seed(1234)

mydata.amount.1 = mydata.amount %>%
  mutate(Exposure.1 = mydata.amount$ClaimAmount*mydata.amount$Exposure/(mydata.amount$ClaimAmount%*%mydata.amount$Exposure))
#%>% select(-Exposure)
L.1 = nrow(mydata.amount.1)

d.1 = sort(sample(L.1, floor(L.1/2)))
amount.train.1 = mydata.amount.1[d.1,]
amount.test.1 = mydata.amount.1[-d.1,]
#glimpse(amount.train.1)
amount.rf = randomForest(ClaimAmount ~. , data = amount.train.1, ntree = 1000,nodesize = 500)

Prédiction sur l’échantillon test

amount.rf.test = predict(amount.rf, newdata = amount.test.1)
summary(amount.rf.test)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1147    1354    1574    2243    1990  513900

Erreur de prédiction

mse.amount.rf = mean((amount.rf.test-amount.test.1$ClaimAmount)^2)
mse.amount.rf
## [1] 120602938

IV. Modélisation par les réseaux de neurones

IV.1 Modélisation de la fréquence

Apprentissage

Prédiction sur l’échantillon de test

freq.nnet.test = predict(freq.nnet, freq.test.1)
summary(freq.nnet.test)
##        V1          
##  Min.   :0.007229  
##  1st Qu.:0.030274  
##  Median :0.040704  
##  Mean   :0.044742  
##  3rd Qu.:0.054472  
##  Max.   :0.317264
#print(freq.nnet.test)
#for(i in 1:length(freq.nnet.test)){if(freq.nnet.test[i] <= 0){freq.nnet.test[i]=0}}

Erreur de prédiction

mse.freq.nnet = mean((freq.nnet.test - freq.test.1$ClaimNb )^2)
mse.freq.nnet
## [1] 0.04052797

IV.2 Modélisation du montant des sinistres

Apprentissage

Prédiction sur l’échantillon de test

amount.nnet.test =max(amount.train$ClaimAmount)* predict(amount.nnet,amount.test )
summary(amount.nnet.test)
##        V1        
##  Min.   :-111.1  
##  1st Qu.:1263.8  
##  Median :1433.7  
##  Mean   :1463.2  
##  3rd Qu.:1642.4  
##  Max.   :3331.2
#max(amount.train$ClaimAmount)*

Erreur de prédiction

mse.amount.nnet = mean((amount.nnet.test-amount.test$ClaimAmount)^2)
mse.amount.nnet
## [1] 3116405

V. Modélisation par le Gradient Boosting Machine

V.1 Modélisation de la fréquence

Nous utiliserons la fonction gbm du package éponyme. Apprentissage

set.seed(1234)

e1 = sort(sample(N, floor(N/5)))
freq.train.1 = mydata.freq[e1,]
freq.test.1 = mydata.freq[-e1,]


freq.gbm = gbm(ClaimNb~., data = freq.train.1,
               distribution = "poisson",n.trees = 2000, 
               shrinkage = 0.1,n.minobsinnode = 1000)
summary(freq.gbm)

##                 var    rel.inf
## Power         Power 20.0712666
## Exposure   Exposure 17.7463596
## DriverAge DriverAge 17.1290939
## Density     Density 15.7034882
## Region       Region 15.1152414
## Brand         Brand  9.2556291
## CarAge       CarAge  4.5049579
## Gas             Gas  0.4739633

Nous utilisons la fonction gbm.perf() pour estimer un nombre d’itérations optimal.

print(best.iter)
## [1] 54

Nous utiliserons donc 54, 54, 54, 54, 54, 54 itérations pour nos prédictions. Prédiction sur l’échantillon de test

freq.gbm.test = predict(freq.gbm,freq.test, n.trees = best.iter)
summary(exp(freq.gbm.test))
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.008284 0.024640 0.042170 0.038910 0.050520 0.128100

Erreur de prédiction

lse.freq.gbm = mean((exp(freq.gbm.test)- freq.test$ClaimNb)^2)
lse.freq.gbm
## [1] 0.04092621

V.2 Modélisation du montant des sinistres

Apprentissage

#amount.train$ClaimAmount = amount.train$ClaimAmount
amount.gbm = gbm(ClaimAmount~., data = amount.train,
               distribution = "gaussian",n.trees = 2000, 
               shrinkage = 0.1,n.minobsinnode = 500, weights = Exposure)
summary(amount.gbm)

##                 var    rel.inf
## Power         Power 31.0463764
## Region       Region 25.9018404
## Brand         Brand 10.9749917
## Density     Density 10.6672954
## Exposure   Exposure  7.4978601
## DriverAge DriverAge  7.4162482
## CarAge       CarAge  5.9020136
## Gas             Gas  0.5933741
best.iter.amount = gbm.perf(amount.gbm,plot.it = TRUE,oobag.curve = TRUE,overlay = TRUE, method="OOB")
## Warning in gbm.perf(amount.gbm, plot.it = TRUE, oobag.curve = TRUE, overlay
## = TRUE, : OOB generally underestimates the optimal number of iterations
## although predictive performance is reasonably competitive. Using cv.folds>0
## when calling gbm usually results in improved predictive performance.

#best.iter = gbm.perf(freq.gbm,plot.it = TRUE,oobag.curve = TRUE,overlay = TRUE,method="cv")
summary(best.iter.amount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      29      29      29      29      29      29

Prédiction sur l’échantillon de test

amount.gbm.test = predict(amount.gbm,amount.test, n.trees = best.iter.amount)
summary(amount.gbm.test)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1156    1319    1425    1429    1535    1827

Erreur de prédiction

lse.amount.gbm = mean((amount.gbm.test- amount.test$ClaimAmount)^2)
lse.amount.gbm
## [1] 3099919

VI. Synthèse des résultats

Nous avons donc mis en oeuvre 4 algorithmes différents de Machine learning. Nous allons maintenant comparer leurs qualités de prédiction et en choisir un seul.