Construir la base de una tarifa de seguros de automóviles

Introducción

Paquetes

library(kableExtra)
library(tidyverse)
library(vcd)
library(MASS)
library(pscl)
library(glmnet)
library(mpath)

La base de datos

Supongamos que disponemos de una BBDD de una compañía de seguros.

La base corresponde a una xtracción aleatoria de 28360 pólizas de la BBDD de experiencia de siniestralidad de una compañía de seguros y consta de 19 variables capturando distintas características de las pólizas:

dim(datos)

## [1] 28360    19

str(datos)

## spc_tbl_ [28,360 × 19] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ npoliza     : num [1:28360] 1202726 877706 446080 106795 735046 ...
##  $ iniciop     : num [1:28360] 2e+07 2e+07 2e+07 2e+07 2e+07 ...
##  $ finalp      : num [1:28360] 2e+07 2e+07 2e+07 2e+07 2e+07 ...
##  $ potencia    : num [1:28360] 218 55 90 120 109 200 75 100 185 50 ...
##  $ valorvehi   : num [1:28360] 7950 1750 1850 2500 1950 ...
##  $ antivehi    : num [1:28360] 4 9 4 10 6 3 1 7 0 17 ...
##  $ tipovehi    : num [1:28360] 100 100 100 100 100 100 100 100 100 100 ...
##  $ plazasvehi  : num [1:28360] 5 5 5 5 5 5 5 5 5 5 ...
##  $ edad        : num [1:28360] 39 62 25 57 50 43 65 25 60 35 ...
##  $ anticarnet  : num [1:28360] 16 29 7 36 23 25 39 7 37 7 ...
##  $ sexocondu   : chr [1:28360] "V" "V" "V" "V" ...
##  $ anticia     : num [1:28360] 0 1 4 11 1 5 3 5 0 0 ...
##  $ codigopostal: num [1:28360] 28003 3202 43500 28730 8206 ...
##  $ usovehi     : num [1:28360] 11000 11000 11000 11000 11000 11000 11000 11000 11000 11000 ...
##  $ cuantiasi   : num [1:28360] 274051 0 247 0 -572 ...
##  $ numerosi    : num [1:28360] 4 0 1 0 1 0 0 0 0 0 ...
##  $ siniestro   : num [1:28360] 1 0 1 0 1 0 0 0 0 0 ...
##  $ prov        : num [1:28360] 28 3 43 28 8 43 25 32 1 7 ...
##  $ ccaa        : num [1:28360] 15 3 8 15 8 8 8 9 1 7 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   npoliza = col_double(),
##   ..   iniciop = col_double(),
##   ..   finalp = col_double(),
##   ..   potencia = col_double(),
##   ..   valorvehi = col_double(),
##   ..   antivehi = col_double(),
##   ..   tipovehi = col_double(),
##   ..   plazasvehi = col_double(),
##   ..   edad = col_double(),
##   ..   anticarnet = col_double(),
##   ..   sexocondu = col_character(),
##   ..   anticia = col_double(),
##   ..   codigopostal = col_double(),
##   ..   usovehi = col_double(),
##   ..   cuantiasi = col_double(),
##   ..   numerosi = col_double(),
##   ..   siniestro = col_double(),
##   ..   prov = col_double(),
##   ..   ccaa = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

npoliza: Número de poliza
iniciop: Inicio del periodo exposición.
finalp: Final del period exposición.
potencia: Potencia del vehículo (en CV).
valorvehi: valor del vehículo (en euros).
antivehi: antigüedad del vehículo (en años).
tipovehi: Tipo de vehículo.
plazasvehi: Plazas del vehículo.
edad: Edad del conductor.
anticarnet: Antigüedad del carnet.
sexocondu: Sexo del conductor.
anticia: Antigüedad en la compañía.
codigopostal: Código postal.
usovehi: Uso del vehículo.
cuantiasi: Pérdidas asociadas a siniestros.
numerosi: Número de siniestros
siniestro: Indicador de si ha declarado o no siniestro
prov: Provincia.
ccaa: Comunidad Autónoma.

¿Cual sería nuestra variable respuesta para nuestro modelo? ¿Qué variable(s) no podemos utilizar como predictor(es)? ¿Por qué? ¿Tenemos más variables de las que aparentemente tenemos?

Tratamiento previo de la base

Muchas de las variables no tienen la clase que les correspondería.

# Descartamos el número de póliza
datos <- subset(datos, select = -c(npoliza))
head(datos)

## # A tibble: 6 × 18
##    iniciop   finalp potencia valorvehi antivehi tipovehi plazasvehi  edad
##      <dbl>    <dbl>    <dbl>     <dbl>    <dbl>    <dbl>      <dbl> <dbl>
## 1 20011015 20020101      218      7950        4      100          5    39
## 2 20000101 20001121       55      1750        9      100          5    62
## 3 20010224 20020101       90      1850        4      100          5    25
## 4 20000101 20000704      120      2500       10      100          5    57
## 5 19990513 20000101      109      1950        6      100          5    50
## 6 20010101 20010511      200      4800        3      100          5    43
## # ℹ 10 more variables: anticarnet <dbl>, sexocondu <chr>, anticia <dbl>,
## #   codigopostal <dbl>, usovehi <dbl>, cuantiasi <dbl>, numerosi <dbl>,
## #   siniestro <dbl>, prov <dbl>, ccaa <dbl>

# Convertimos en formato fecha: as.Date(x, "%m/%d/%Y")  
datos$iniciop <- paste0(substr(datos$iniciop, 5, 6), "/",
                        substr(datos$iniciop, 7, 8), "/",
                        substr(datos$iniciop, 1, 4) )
datos$iniciop <- as.Date(datos$iniciop , "%m/%d/%Y")
datos$finalp <- paste0(substr(datos$ finalp, 5, 6), "/",
                       substr(datos$ finalp, 7, 8), "/",
                       substr(datos$ finalp, 1, 4) )
datos$finalp <- as.Date(datos$ finalp , "%m/%d/%Y")

# Observamos cambios  
head(datos)

## # A tibble: 6 × 18
##   iniciop    finalp     potencia valorvehi antivehi tipovehi plazasvehi  edad
##   <date>     <date>        <dbl>     <dbl>    <dbl>    <dbl>      <dbl> <dbl>
## 1 2001-10-15 2002-01-01      218      7950        4      100          5    39
## 2 2000-01-01 2000-11-21       55      1750        9      100          5    62
## 3 2001-02-24 2002-01-01       90      1850        4      100          5    25
## 4 2000-01-01 2000-07-04      120      2500       10      100          5    57
## 5 1999-05-13 2000-01-01      109      1950        6      100          5    50
## 6 2001-01-01 2001-05-11      200      4800        3      100          5    43
## # ℹ 10 more variables: anticarnet <dbl>, sexocondu <chr>, anticia <dbl>,
## #   codigopostal <dbl>, usovehi <dbl>, cuantiasi <dbl>, numerosi <dbl>,
## #   siniestro <dbl>, prov <dbl>, ccaa <dbl>

str(datos)

## tibble [28,360 × 18] (S3: tbl_df/tbl/data.frame)
##  $ iniciop     : Date[1:28360], format: "2001-10-15" "2000-01-01" ...
##  $ finalp      : Date[1:28360], format: "2002-01-01" "2000-11-21" ...
##  $ potencia    : num [1:28360] 218 55 90 120 109 200 75 100 185 50 ...
##  $ valorvehi   : num [1:28360] 7950 1750 1850 2500 1950 ...
##  $ antivehi    : num [1:28360] 4 9 4 10 6 3 1 7 0 17 ...
##  $ tipovehi    : num [1:28360] 100 100 100 100 100 100 100 100 100 100 ...
##  $ plazasvehi  : num [1:28360] 5 5 5 5 5 5 5 5 5 5 ...
##  $ edad        : num [1:28360] 39 62 25 57 50 43 65 25 60 35 ...
##  $ anticarnet  : num [1:28360] 16 29 7 36 23 25 39 7 37 7 ...
##  $ sexocondu   : chr [1:28360] "V" "V" "V" "V" ...
##  $ anticia     : num [1:28360] 0 1 4 11 1 5 3 5 0 0 ...
##  $ codigopostal: num [1:28360] 28003 3202 43500 28730 8206 ...
##  $ usovehi     : num [1:28360] 11000 11000 11000 11000 11000 11000 11000 11000 11000 11000 ...
##  $ cuantiasi   : num [1:28360] 274051 0 247 0 -572 ...
##  $ numerosi    : num [1:28360] 4 0 1 0 1 0 0 0 0 0 ...
##  $ siniestro   : num [1:28360] 1 0 1 0 1 0 0 0 0 0 ...
##  $ prov        : num [1:28360] 28 3 43 28 8 43 25 32 1 7 ...
##  $ ccaa        : num [1:28360] 15 3 8 15 8 8 8 9 1 7 ...

# Calculamos algunos resumenes
apply(datos[,c(3, 5:11, 13, 15:18)], 2, table)

## $potencia
## 
##   17   29   30   31   32   34   35   37   38   39   40   41   42   43   44   45 
##    2   10    5   25   20   42   43    3   36    6  250    6  201  146   45  789 
##   46   47   48   49   50   51   52   53   54   55   56   57   58   59   60   61 
##   65   19   74   80  440   45  147   70  491 1366    2   70  260   14 3019   50 
##   62   63   64   65   66   67   68   69   70   71   72   73   74   75   76   77 
##   17  351  234 1027    4  121  310  112  464  636  395   39   80 1826   91   12 
##   78   79   80   81   82   83   84   85   86   87   88   89   90   91   92   93 
##  126   76  667    1  569  175   79  215  144   43  286   59 2626    4  611   22 
##   94   95   96   97   98   99  100  101  102  103  104  105  106  107  108  109 
##  176  362   45   37   11   87  957   21  141  167   31  298   32  108   56   64 
##  110  111  112  113  114  115  116  117  118  119  120  121  122  123  124  125 
##  616   53  429  119  113 1182   33   66   85    6  275   26  151   64   27  477 
##  126  128  129  130  131  132  133  134  135  136  137  138  139  140  141  142 
##   35   30  106  154   24   25   92    6   34  380   56   80   32  133    6    5 
##  143  144  145  146  147  148  150  151  153  154  155  156  157  158  159  160 
##   82    7   22    5   27    4  401   12    1    2   29    4   10    4    1   58 
##  161  162  163  165  166  167  168  169  170  171  173  174  175  177  178  180 
##    6    3   17   54   10   12    2   10  129    5    1   46   20   26    7   25 
##  181  182  184  185  186  187  188  190  192  193  194  195  197  200  201  203 
##    5    2   32   58    1    3   39    6   40   69    7    2    9   28    1    8 
##  204  205  207  208  210  211  215  216  218  220  221  222  223  224  225  227 
##   21    3    4    2    3   17    9    2   40   29    8    1    3    6   10    4 
##  228  230  231  235  237  238  240  245  248  250  252  254  272  286  291  295 
##    4    9    9    3    2   17    6    7    1    2    1    3    1    1    1    2 
##  300 
##    1 
## 
## $antivehi
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
##  969 1996 2233 2277 1847 1566 1471 1614 1601 1828 1884 1807 1728 1369 1011  659 
##   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   32 
##  436  449  352  392  153  141  119   82   86   61   61   62   52   43    9    2 
## 
## $tipovehi
## 
##   100   120   150   200   250   300   310 
## 25019   349  1281  1199   137    50   317 
## 
## $plazasvehi
## 
##     5     6 
## 28097   263 
## 
## $edad
## 
##  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37 
##   6  12  16  36  70 105 108 185 263 339 350 401 454 562 618 691 728 735 685 654 
##  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57 
## 753 731 815 800 763 771 777 836 816 890 874 851 899 865 843 844 786 738 646 637 
##  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77 
## 531 471 477 406 372 357 332 339 247 249 229 212 166 135 140 107 104  94  89  74 
##  78  79  80  81  82  83  84  85  86  87  88  89  90 
##  70  49  45  31  25  15  12   7   7   2   6   3   4 
## 
## $anticarnet
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
##   38   44   68  118  190  229  277  413  527  634  692  742  836  829  868  843 
##   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31 
##  817  791  879  894 1015 1140 1150 1083 1106 1029 1047  999  955  880  868  862 
##   32   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47 
##  741  762  685  604  573  424  378  266  212  183  145   81   72   71   45   41 
##   48   49   50   51   52   53   54   55   56   57   58   59   60   61   62   63 
##   40   23   26   23   15   19   11   13   11    4    4    4    5    3    4    3 
##   64   65   68   71   72 
##    1    2    1    1    1 
## 
## $sexocondu
## 
##     M     V 
##  6008 22352 
## 
## $anticia
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 6008 5079 4172 3414 2080 1280  907  771  685  601  513  502  479  469  355  322 
##   16   17   18 
##  277  247  199 
## 
## $usovehi
## 
##     0 11000 11100 13100 13300 13500 14100 14200 15000 21000 21001 21002 21003 
##    11 25852   869    55    32     9    10     1    11    59   581     4   490 
## 21004 21005 21007 21102 23100 23201 23202 
##    23   179   169     1     1     2     1 
## 
## $numerosi
## 
##     0     1     2     3     4     5     6     7     8     9 
## 20447  5463  1408   591   282   108    40    13     5     3 
## 
## $siniestro
## 
##     0     1 
## 20447  7913 
## 
## $prov
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
##  446  151 1524  298   23  302 1136 5961  369  157  439 1033  135  370  920   93 
##   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 
##  571  590   36  318  100  184  420  217  777  201  124  832 1111  549  279  242 
##   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47   48 
##  700   88  507  652  203  968  327   83  549  160  495  112  348 1884  301  358 
##   49   50 
##   88  629 
## 
## $ccaa
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 1122  763 4441 3877 1532  459 1136 7804 1938  925  549  279  327  700  832  201 
##   17 
## 1475

# Declaramos como factores las variables que corresponda  
datos$tipovehi <- as.factor(datos$tipovehi)
datos$codigopostal <- as.factor(datos$codigopostal)
datos$prov <- as.factor(datos$prov)
datos$ccaa <- as.factor(datos$ccaa)
str(datos)

## tibble [28,360 × 18] (S3: tbl_df/tbl/data.frame)
##  $ iniciop     : Date[1:28360], format: "2001-10-15" "2000-01-01" ...
##  $ finalp      : Date[1:28360], format: "2002-01-01" "2000-11-21" ...
##  $ potencia    : num [1:28360] 218 55 90 120 109 200 75 100 185 50 ...
##  $ valorvehi   : num [1:28360] 7950 1750 1850 2500 1950 ...
##  $ antivehi    : num [1:28360] 4 9 4 10 6 3 1 7 0 17 ...
##  $ tipovehi    : Factor w/ 7 levels "100","120","150",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ plazasvehi  : num [1:28360] 5 5 5 5 5 5 5 5 5 5 ...
##  $ edad        : num [1:28360] 39 62 25 57 50 43 65 25 60 35 ...
##  $ anticarnet  : num [1:28360] 16 29 7 36 23 25 39 7 37 7 ...
##  $ sexocondu   : chr [1:28360] "V" "V" "V" "V" ...
##  $ anticia     : num [1:28360] 0 1 4 11 1 5 3 5 0 0 ...
##  $ codigopostal: Factor w/ 4422 levels "1001","1002",..: 2242 141 3767 2340 589 3729 2064 2768 1 413 ...
##  $ usovehi     : num [1:28360] 11000 11000 11000 11000 11000 11000 11000 11000 11000 11000 ...
##  $ cuantiasi   : num [1:28360] 274051 0 247 0 -572 ...
##  $ numerosi    : num [1:28360] 4 0 1 0 1 0 0 0 0 0 ...
##  $ siniestro   : num [1:28360] 1 0 1 0 1 0 0 0 0 0 ...
##  $ prov        : Factor w/ 50 levels "1","2","3","4",..: 28 3 43 28 8 43 25 32 1 7 ...
##  $ ccaa        : Factor w/ 17 levels "1","2","3","4",..: 15 3 8 15 8 8 8 9 1 7 ...

# Colapsamos usos de vehiculo
datos$usovehi2 <- ifelse(datos$usovehi < 20000, "G1", "G2")

# Descartamos uso de vehiculo, poca varianza
datos <- subset(datos, select = -c(usovehi))

Descartamos siniestro y cuantiasi como predictores, ¿por qué?

datos <- subset(datos, select = -c(cuantiasi, siniestro))

Bases de Training y de Test

Recuerda: No se deben "mirar" los datos de test para especificar la estrategia/el modelo de predicción, al fin y al cabo estos no se conocen en la realidad.

Partimos al azar base en dos: 66.67% entrenamiento (18.907 casos) y 33.33% test (9.453 casos).

# Permutamos los datos al azar
set.seed(216514)
azar <- sample(1:nrow(datos))
casos <- round(nrow(datos)*0.66667)
train <- datos[azar[1:casos],]
test <- datos[azar[(casos+1):nrow(datos)],]
Tabla.Training <- table(train$numerosi)
Tabla.Test <- table(test$numerosi)

Tabla.Training %>%
  kbl(col.names = c("Número Reclamaciones", "Frecuencia"),
      align = rep('c', 2)) %>%
  kable_styling()

Número Reclamaciones	Frecuencia
0	13531
1	3706
2	949
3	419
4	192
5	72
6	24
7	9
8	3
9	2

Tabla.Test %>%
  kbl(col.names = c("Número Reclamaciones", "Frecuencia"),
      align = rep('c', 2)) %>%
  kable_styling()

Número Reclamaciones	Frecuencia
0	6916
1	1757
2	459
3	172
4	90
5	36
6	16
7	4
8	2
9	1

¿Qué alternativa de partición de la base propondrías?

Análisis Exploratorio de Datos

Variable respuesta

Comenzamos a conocer la variable respuesta

barplot(table(train$numerosi))

¿Sobredispersión? ¿Exceso de ceros?

summary(train$numerosi)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.4357  1.0000  9.0000

var(train$numerosi)

## [1] 0.751857

Comparación con una Poisson

observado <- table(train$numerosi)
predict0 <- dpois(0:9, mean(train$numerosi))*nrow(train)
barplot(rbind(observado, predict0), las = 1,
        beside = TRUE, names.arg = c(0:8, ">=9"),
        col=c("blue", "green"))

Comparación con una Binomial Negativa

fit.bn <- goodfit(train$numerosi, type = "nbinomial", 
                  method = "MinChisq")
barplot(rbind(observado, fit.bn$fitted[1:10]), las = 1,
        beside = TRUE, names.arg = c(0:8, ">9"),
        col=c("blue", "green"))

Relaciones variable respuesta-predictores

Tiempo de exposición al riesgo

¿Cómo debemos tratarlo? ¿cómo un predictor? ¿cómo un peso? ¿como offset?

# Tiempo de exposición (Predictor, Peso u Offset)
train$tiempo <- as.numeric(difftime(train$finalp,                
                                    train$iniciop, units="days"))
summary(train$tiempo)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    88.0   185.0   183.6   277.0   366.0

¿Qué hacemos con los que están cero días expuestos al riesgo? ¿Cómo es posible?

sum(train$tiempo==0)

## [1] 360

plot(density(train$tiempo))

¿Cómo son esas observaciones?

table(train$numerosi[train$tiempo==0])

## 
##   0   1   5 
## 355   4   1

train[train$tiempo==0 & train$numerosi==1,]

## # A tibble: 4 × 17
##   iniciop    finalp     potencia valorvehi antivehi tipovehi plazasvehi  edad
##   <date>     <date>        <dbl>     <dbl>    <dbl> <fct>         <dbl> <dbl>
## 1 1999-09-30 1999-09-30       65      2050        1 100               5    27
## 2 2000-04-19 2000-04-19      135      3750        0 100               5    34
## 3 1999-04-06 1999-04-06       72      2250       14 310               5    42
## 4 1999-01-19 1999-01-19       55      1350       14 100               5    34
## # ℹ 9 more variables: anticarnet <dbl>, sexocondu <chr>, anticia <dbl>,
## #   codigopostal <fct>, numerosi <dbl>, prov <fct>, ccaa <fct>, usovehi2 <chr>,
## #   tiempo <dbl>

train <- train[train$tiempo != 0, ]
cor(train$tiempo, train$numerosi)

## [1] 0.2804692

plot(sort(unique(train$tiempo)), 
     tapply(train$numerosi, train$tiempo, mean))

Potencia del vehículo

summary(train$potencia)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   29.00   60.00   80.00   86.59  105.00  300.00

plot(density(train$potencia, bw=8))

medias <- tapply(train$numerosi, train$potencia, mean)
plot(sort(unique(train$potencia)), medias, pch=20)

cor(train$potencia, train$numerosi)

## [1] 0.05695214

¿Existe relación con la variable respuesta? ¿Es lineal?

¿Y si la discretizamos? ¿Cómo? ¿Cuántos grupos? ¿Dónde cortes?

Potencia discretizada

# Potencia Discretizada 
qnt <- quantile(train$potencia,seq(0,1,.2))
qnt[1] <- qnt[1]-0.01; qnt[5] <- qnt[5]+0.01
train$g.potencia <- cut(train$potencia,  breaks= qnt)
medias <- tapply(train$numerosi, train$g.potencia, mean)
plot(medias, pch=20, xaxt= "n", xlab = "") # log(medias)
axis(1, at=1:5, labels=levels(train$g.potencia))

¿Construimos otros predictores basados en fechas?

% de tiempo expuesto en cada mes
% de tiempo expuesto en cada estación
% de tiempo expuesto de Navidad, en Semana Santa,…

¿Discretizamos la variable?

Valor del vehículo

plot(density(train$valorvehi))

plot(sort(unique(train$valorvehi)), 
     tapply(train$numerosi, train$valorvehi, mean))

¿Y si discretizamos?

Discretizar parace que permite ver las relaciones de una manera más nítida, es algo que vamos a hacer con mucha frecuencia.… construyamos una función.….

# Función pintar y discretizar
plot.dist <- function(predictor, cortes=8, variable=F){
      qnt <- quantile(predictor, seq(0,1,1/cortes))
      qnt[1] <- qnt[1]-0.01; qnt[cortes] <- qnt[cortes]+0.01
      V.discreta <- cut (predictor,  breaks= qnt)
      medias <- tapply(train$numerosi, V.discreta, mean)
      plot(medias, pch=20, xaxt= "n", xlab = "")
      axis(1, at = 1:cortes, labels=levels(V.discreta))
      if (variable == T) return(V.discreta)
} # Fin función

Valor del vehículo discretizada

plot.dist(train$valorvehi)

train$g.valorvehi <- plot.dist(train$valorvehi, 5, TRUE)

cor(train$valorvehi, train$numerosi)

## [1] 0.06386046

Antigüedad del vehículo

plot(density(train$antivehi, bw = .5))

plot(sort(unique(train$antivehi)), 
     tapply(train$numerosi, train$antivehi, mean))

plot.dist(train$antivehi, 9)

cor(train$antivehi, train$numerosi)

## [1] -0.1280255

Tipo de vehículo

Esta variable es de tipo factor, es preciso otro tipo de análisis.

summary(train$tipovehi)

##   100   120   150   200   250   300   310  NA's 
## 16397   230   865   730    82    26   211     6

plot(levels(train$tipovehi), 
     tapply(train$numerosi, train$tipovehi, mean))

mod1 <- lm(train$numerosi ~ train$tipovehi)
anova(mod1)

## Analysis of Variance Table
## 
## Response: train$numerosi
##                   Df  Sum Sq Mean Sq F value Pr(>F)   
## train$tipovehi     6    15.3 2.54615  3.3451 0.0027 **
## Residuals      18534 14107.4 0.76116                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

confint(mod1)

##                         2.5 %      97.5 %
## (Intercept)        0.43166578  0.45837508
## train$tipovehi120  0.07186739  0.29896132
## train$tipovehi150 -0.05959191  0.05972445
## train$tipovehi200 -0.15765185 -0.02827942
## train$tipovehi250 -0.26848438  0.11015083
## train$tipovehi300 -0.20373609  0.46754138
## train$tipovehi310 -0.09430745  0.14265522

Plazas de vehículo

No presenta mucha variabilidad, ¿qué tal su impacto?

```{rtable(train$plazasvehi)}



```r
tapply(train$numerosi, train$plazasvehi, mean)

##         5         6 
## 0.4442329 0.3832335

Edad

plot(density(train$edad, bw=.5))

plot(sort(unique(train$edad)), 
     tapply(train$numerosi, train$edad, mean))

plot.dist(train$edad, 9)

cor(train$edad, train$numerosi)

## [1] -0.04166107

Antigüedad del carnet

plot(density(train$anticarnet, bw=.5))

plot(sort(unique(train$anticarnet)), 
     tapply(train$numerosi, train$anticarnet, mean))

plot.dist(train$anticarnet, 12)

cor(train$anticarnet, train$numerosi)

## [1] -0.03409736

Sexo conductor

table(train$sexocondu)

## 
##     M     V 
##  3904 14643

tapply(train$numerosi, train$sexocondu, mean)

##         M         V 
## 0.5804303 0.4072253

Antigüedad en la compañía

plot(density(train$anticia, bw = .5))

plot(sort(unique(train$anticia)), 
     tapply(train$numerosi, train$anticia, mean))

plot.dist(train$anticia, 4)

cor(train$anticia, train$numerosi)

## [1] -0.08779265

train$g.anticia <- ifelse(train$anticia < 3, "New", "Vet")

Usos de Vehículo

table(train$usovehi2)

## 
##    G1    G2 
## 17616   931

tapply(train$numerosi, train$usovehi2, mean)

##        G1        G2 
## 0.4472071 0.3770140

Algunas conclusiones del EDA

Casi todas las variables estudiadas y las transformaciones construidas muestran (individualmente) capacidad predictiva (salvo plazas del vehículo y tipo de vehículo). Se incluirán inicialmente todas ellas. En el caso de la edad (y otras, e.g. antigüedad del carnet), hemos observado relación no lineal.

Incluimos edad^2 y edad^3, también como predictores (alternativa modelos GAM).

# Edad cuadrado y cubo
train$edad2 <- train$edad^2
train$edad3 <- train$edad^3

¿Más predictores?

Obviamente es posible construir más predictores, ideas:

Discretizando más variables y/o de otro modo.
Incorporando efectos de interacción.
Construyendo nuevas variables a partir de las ya disponibles.
Derivando nuevas variables a partir de las ya disponibles: por ejemplo, a partir del código postal.
...

Propuestas de modelización

Especificación:

Variables que van a entrar en el predictor lineal.
Distribución para la variable respuesta
Función link
Y con enfoque ML, tipo de regularización/penalización

Distribuciones para la variable respuesta

Como modelos para la variable respuesta podemos testar:

Poisson
Binomial Negativa
ZI (inflada en cero) + Poisson
ZI (inflada en cero) + BN
Normal

Funciones link

Como funciones link:

Logaritmo
Identidad

Regularización

Penalización: * Net Elastic (Lasso, Ridge)

Modelos sin regularización

Cuando utilizamos modelos sin regularización hemos de estar atentos a los problemas de multicolinealidad

Antes de comenzar a modelizar hemos de eliminar algunas variables por que es lo razonable, lógico, desde el punto de vista de modelización y de computación.

Eliminamos codigopostal, ccaa, uso de vehiculo, fechas inicio y final ¿Por qué?

train <- subset(train, select = -c(codigopostal, ccaa,
                                   iniciop, finalp))

Ajuste modelos (clásicos) de Poisson

Observa la sintaxis

Po.log <- glm(numerosi ~ . , data = train, 
              family=poisson(link="log"))
starting.values <- coef(Po.log)
starting.values[starting.values < 0] <- 0.0001
Po.id <- glm(numerosi ~ . , data= train,
             start=starting.values,
             family=poisson(link="identity"))

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning in log(y/mu): Se han producido NaNs

## Warning: step size truncated due to divergence

## Warning: step size truncated: out of bounds

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm stopped at boundary value

Ajuste modelo (clásicos) Binomial Negativo

# library(MASS)
BN.log <- glm.nb(numerosi ~ . , data= train, link = log)
# starting.values <- coef(BN.log)
# starting.values[starting.values<0] <- 0.000001
# BN.id <- glm.nb(numerosi ~ . , data= train,
#                start= starting.values, link = identity)

Ajustes modelos (clásicos) inflados en cero

# library(pscl)
ZI.Po <- zeroinfl(numerosi ~ ., data=train, dist="poisson")
ZI.BN <- zeroinfl(numerosi ~ ., data=train, dist="negbin")

NORMAL

No.id <- glm(numerosi ~ ., data=train)
No.log <- glm(numerosi ~ ., data=train, start = coef(No.id),
              family=gaussian(link="log"))

¿Tiene sentido ajustar con respuesta Normal?

El rol del tiempo de exposición: Modelos con Offsets (y pesos)

Poisson con offset

Po.log.off <- glm(numerosi ~ . ,
                data = subset(train , select = -c(tiempo)),
                offset = log(train$tiempo/365.25),  
                family = poisson(link="log"))

Poisson con pesos

Po.log.wgt <- glm(numerosi ~ . ,
                data = subset(train , select = -c(tiempo)),
                weight = train$tiempo/365.25,  
                family = poisson(link="log"))

Comparamos AIC's

AIC(Po.log.off)

## [1] 31281.97

AIC(Po.log.wgt)*nrow(train)/sum(train$tiempo/365.25)

## [1] 39318.87

Poisson con offset (especificacion alternativa)

of_var <- log(train$tiempo/365.25)
formula <- numerosi ~ potencia + valorvehi + antivehi +
             tipovehi + plazasvehi + edad + anticarnet +
             sexocondu + anticia + prov + usovehi2 +
             g.potencia + offset(of_var) +
             g.valorvehi + g.anticia + edad2 + edad3

Po.log.off2 <- glm(formula, data = train, family = poisson)

Binomial Negativa con offset

BN.log.off <- glm.nb(formula, data = train, link = log)

Modelos inflados de ceros, con offset

formula2 <- formula(numerosi ~ potencia + valorvehi +
             antivehi + tipovehi + plazasvehi + edad + 
             anticarnet + sexocondu + anticia + 
             prov + usovehi2 + g.potencia +
             offset(log(tiempo/365.25)) +
             g.valorvehi + g.anticia + edad2 + edad3
           | potencia + valorvehi + antivehi + tipovehi +
             plazasvehi + edad + anticarnet + sexocondu +
             anticia + prov + usovehi2 + g.potencia +
             g.valorvehi + g.anticia + edad2 + edad3)

ZI.Po.off <- zeroinfl(formula2,
                      data = train, dist="poisson")

ZI.BN.off <- zeroinfl(formula2,
                      data = train, dist="negbin")

Evaluando modelos. Evaluación clásica

# EVALUACION CLASICA DE MODELOS
AICs <- c(AIC(Po.log), AIC(Po.id), AIC(BN.log), AIC(ZI.Po), AIC(ZI.BN),
          AIC(No.id), AIC(No.log), AIC(Po.log.off), 
          # ‘Ajuste’ por número de observaciones
          AIC(Po.log.wgt)*nrow(train)/sum(train$tiempo/365.25),
          AIC(BN.log.off), AIC(ZI.Po.off), AIC(ZI.BN.off))

errores_clasicos <- data.frame(
  modelo = c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", 
             "No.id", "No.log", "Po.log.off", "Po.log.wgt", 
             "BN.log.off", "ZI.Po.off", "ZI.BN.off"),
  AICs)

errores_clasicos %>%
  kbl(align = rep('c', 3)) %>%
  kable_styling()

modelo	AICs
Po.log	31533.58
Po.id	49633.67
BN.log	30508.43
ZI.Po	30520.97
ZI.BN	30186.29
No.id	45478.90
No.log	45370.31
Po.log.off	31281.97
Po.log.wgt	39318.87
BN.log.off	30297.79
ZI.Po.off	30749.12
ZI.BN.off	30230.62

¿Qué modelo elegirías? ¿Por qué?

Evaluando modelos. Evaluación enfoque ML

Definicion nuevos predictores definidos en el conjunto de training en test

qnt <- quantile(train$potencia,seq(0,1,.2))
qnt[1] <- qnt[1]-0.01; qnt[5] <- qnt[5]+0.01
test$g.potencia <- cut(test$potencia,  breaks= qnt)
test$tiempo <- as.numeric(difftime(test$finalp,                
                                   test$iniciop,units="days"))

dist.test <- function(predictor.train, predictor.test, cortes=10){
       qnt <- quantile(predictor.train, seq(0,1,1/cortes))
       qnt[1] <- qnt[1]-0.01; qnt[cortes] <- qnt[cortes]+0.01
       V.discreta <- cut (predictor.test,  breaks= qnt)
       return(V.discreta)
}

test$g.valorvehi <- dist.test(train$valorvehi, test$valorvehi, 5)
test$g.anticia <- ifelse(test$anticia < 3, "New", "Vet")
test$edad2 <- test$edad^2
test$edad3 <- test$edad^3
test <- subset(test, select = -c(codigopostal, ccaa, iniciop, finalp))
test <- test[test$tiempo != 0,]

Eliminamos del conjunto test los casos, caso de existir, para los que no se puede crear un predictor

test0 <- test[rowSums(is.na(test))==0,]

Predicciones en el conjunto de test (modelos sin regularización)

Observa la sintaxis para obtener las predicciones.

p.Po.log <- predict(Po.log, test0 , type="response")
p.Po.id <- predict(Po.id, test0 , type="response")
p.BN.log <- predict(BN.log, test0 , type="response")
p.ZI.Po <- predict(ZI.Po, test0 , type="response")
p.ZI.BN <- predict(ZI.BN, test0 , type="response")
p.No.id <- predict(No.id, test0)
p.No.log <- predict(No.log, test0 , type="response")
nd <- data.frame(test0, of_var = log(test0$tiempo/365.25)) 
p.Po.log.off <- predict(Po.log.off2, nd , type = "response")
p.BN.log.off <- predict(BN.log.off, nd, type = "response")
p.ZI.Po.off <- predict(ZI.Po.off, test0 , type="response")
p.ZI.BN.off <- predict(ZI.BN.off, test0 , type="response")
# p.Po.log.wgt <- predict(Po.log.wgt, test0 , type="response")}

Errores cuadráticos

# Función de perdida cuadrática
e.Po.log <- sum((p.Po.log - test0$numerosi)^2)
e.Po.id <- sum((p.Po.id - test0$numerosi)^2)
e.BN.log <- sum((p.BN.log - test0$numerosi)^2)
e.ZI.Po <- sum((p.ZI.Po - test0$numerosi)^2)
e.ZI.BN <- sum((p.ZI.BN - test0$numerosi)^2)
e.No.id <- sum((p.No.id - test0$numerosi)^2)
e.No.log <- sum((p.No.log - test0$numerosi)^2)
e.Po.log.off <- sum((p.Po.log.off - test0$numerosi)^2)
e.BN.log.off <- sum((p.BN.log.off - test0$numerosi)^2)
e.ZI.Po.off <- sum((p.ZI.Po.off - test0$numerosi)^2)
e.ZI.BN.off <- sum((p.ZI.BN.off - test0$numerosi)^2)

e.base <- sum((mean(train$numerosi) - test0$numerosi)^2)

errores_cuadraticos <- c(e.Po.log, e.Po.id, e.BN.log, e.ZI.Po, e.ZI.BN,
                         e.No.id, e.No.log, e.Po.log.off, e.BN.log.off,
                         e.ZI.Po.off, e.ZI.BN.off, e.base)

Errores valor absoluto

a.Po.log <- sum(abs(p.Po.log - test0$numerosi))
a.Po.id <- sum(abs(p.Po.id - test0$numerosi))
a.BN.log <- sum(abs(p.BN.log - test0$numerosi))
a.ZI.Po <- sum(abs(p.ZI.Po - test0$numerosi))
a.ZI.BN <- sum(abs(p.ZI.BN - test0$numerosi))
a.No.id <- sum(abs(p.No.id - test0$numerosi))
a.No.log <- sum(abs(p.No.log - test0$numerosi))
a.Po.log.off <- sum(abs(p.Po.log.off - test0$numerosi))
a.BN.log.off <- sum(abs(p.BN.log.off - test0$numerosi))
a.ZI.Po.off <- sum(abs(p.ZI.Po.off - test0$numerosi))
a.ZI.BN.off <- sum(abs(p.ZI.BN.off - test0$numerosi))

a.base <- sum(abs(mean(train$numerosi) - test0$numerosi))

errores_abs <- c(a.Po.log, a.Po.id, a.BN.log, a.ZI.Po, a.ZI.BN, 
                 a.No.id, a.No.log, a.Po.log.off, a.BN.log.off, 
                 a.ZI.Po.off, a.ZI.BN.off, a.base)

Resumen errores (modelos sin regularización)

errores_sin <- data.frame(
  modelo = c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", 
             "No.id", "No.log", "Po.log.off", "BN.log.off", 
             "ZI.Po.off", "ZI.BN.off", "Null model"),
  errores_cuadraticos,
  errores_abs)

errores_sin[, 2:3] <- round(errores_sin[, 2:3], 2)
errores_sin %>%
  kbl(align = rep('c', 3)) %>%
  kable_styling()

modelo	errores_cuadraticos	errores_abs
Po.log	6006.79	4872.91
Po.id	16994.70	10580.04
BN.log	6023.33	4870.18
ZI.Po	6026.83	4846.27
ZI.BN	6028.87	4845.07
No.id	6078.13	4962.18
No.log	6015.23	4949.90
Po.log.off	6002.92	4825.14
BN.log.off	6010.26	4833.45
ZI.Po.off	6003.26	4815.73
ZI.BN.off	6046.58	4834.45
Null model	6836.01	5728.05

¿Qué modelo elegirías? ¿Por qué? ¿En que se diferencia de una aproximación clásica?

Modelos con regularización

Los Paquetes glmnet y mpath

Cuestiones sobre las funciones

Conversión de data.frame en matriz (construyendo variables dummy para factores e intercepto), necesario para predecir.

model.matrix(~., data)

Modelo con respuesta Poisson

cv.glmnet(x, y, family=“poisson”, alpha, nfolds) glmnet(x, y, family=“poisson”, alpha) cvglmreg(formula, data, family=“poisson”, penalty, nfolds) glmreg(formula, data, family=“poisson”, penalty)

Modelo con respuesta Binomial Negativa

cvglmregNB(formula, data, penalty, nfolds) glmregNB(formula, data, penalty)

Modelos inflados de ceros

cvzipath(formula, data, family, link, penalty, nfolds) zipath(formula, data, family, link, penalty)

Predicción

predict(object, newdata, type) predict(object, newx, which, type)

Más info

https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#poi https://cran.r-project.org/web/packages/mpath/mpath.pdf

Conversión a objetos matriz de predictores

# Convertimos los data.frame en objetos matriz respetando factores
train.m <- model.matrix(~., data=train)
train.m <- train.m[,colnames(train.m) != "numerosi"]

test.m <- model.matrix(~., data=test0)
test.m <- test.m[,colnames(test.m) != "numerosi"]

Modelos con regularización

Se estima un modelo, existe una solucion, por cada valor del parametro de penalización.

Ajustes con glmnet

Observa la sintaxis

Respuesta Normal

#library(glmnet)
#No.re <- glmnet(x = train.m, y = train$numerosi, alpha=0.5)

Respuesta Poisson

#Po.re <- glmnet(x = train.m, y = train$numerosi,
#                family = "poisson", alpha = 0.5)

Respuesta Poisson con offset

#Po.re.off <- glmnet(x = train.m, y = train$numerosi, alpha = 0.5,  
#                   offset = log(train$tiempo/365.25), family = "poisson")

Ajustes con mpath

Observa la sintaxis

Respuesta Poisson

# library(mpath)
#Po.re2 <- glmreg(numerosi ~ . , data = train,
#                 family = "poisson", penalty="enet")

Respuesta Binomial Negativa

#BN.re <- glmregNB(numerosi ~ . , data = train, penalty="enet")

Construir la base de una tarifa de seguros de automóviles

Miguel Ibañez Ureña

2024-06-14

Introducción

Paquetes

La base de datos

Tratamiento previo de la base

Bases de Training y de Test

Análisis Exploratorio de Datos

Variable respuesta

¿Sobredispersión? ¿Exceso de ceros?

Comparación con una Poisson

Comparación con una Binomial Negativa

Relaciones variable respuesta-predictores

Tiempo de exposición al riesgo

Potencia del vehículo

Potencia discretizada

Valor del vehículo

Valor del vehículo discretizada

Antigüedad del vehículo

Tipo de vehículo

Plazas de vehículo

Edad

Antigüedad del carnet

Sexo conductor

Antigüedad en la compañía

Usos de Vehículo

Algunas conclusiones del EDA

¿Más predictores?

Propuestas de modelización

Distribuciones para la variable respuesta

Funciones link

Regularización

Modelos sin regularización

Ajuste modelos (clásicos) de Poisson

Ajuste modelo (clásicos) Binomial Negativo

Ajustes modelos (clásicos) inflados en cero

NORMAL

El rol del tiempo de exposición: Modelos con Offsets (y pesos)

Poisson con offset

Poisson con pesos

Comparamos AIC's

Poisson con offset (especificacion alternativa)

Binomial Negativa con offset

Modelos inflados de ceros, con offset

Evaluando modelos. Evaluación clásica

Evaluando modelos. Evaluación enfoque ML

Definicion nuevos predictores definidos en el conjunto de training en test

Predicciones en el conjunto de test (modelos sin regularización)

Errores cuadráticos

Errores valor absoluto

Resumen errores (modelos sin regularización)

Modelos con regularización

Los Paquetes glmnet y mpath

Cuestiones sobre las funciones

Más info

Conversión a objetos matriz de predictores

Modelos con regularización

Ajustes con glmnet

Respuesta Normal

Respuesta Poisson

Respuesta Poisson con offset

Ajustes con mpath

Respuesta Poisson

Respuesta Binomial Negativa