#Instalar y cargar la librería neuralnet

# Instala e carga la librería neuralnet

#install.packages("neuralnet")
library(neuralnet)
library(readxl)
#install.packages("matrixStats")
library(matrixStats)
#install.packages("mice")
library(mice)
## 
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:matrixStats':
## 
##     count
## The following object is masked from 'package:neuralnet':
## 
##     compute
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#install.packages("recipes")
library(recipes)
## 
## Attaching package: 'recipes'
## The following object is masked from 'package:stats':
## 
##     step
#install.packages("h2o")
library(h2o)
## 
## ----------------------------------------------------------------------
## 
## Your next step is to start H2O:
##     > h2o.init()
## 
## For H2O package documentation, ask for help:
##     > ??h2o
## 
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit https://docs.h2o.ai
## 
## ----------------------------------------------------------------------
## 
## Attaching package: 'h2o'
## The following objects are masked from 'package:stats':
## 
##     cor, sd, var
## The following objects are masked from 'package:base':
## 
##     %*%, %in%, &&, ||, apply, as.factor, as.numeric, colnames,
##     colnames<-, ifelse, is.character, is.factor, is.numeric, log,
##     log10, log1p, log2, round, signif, trunc

Cargar los datos

set.seed(42)
datos <- read_excel("D:/MANADINE/2do año/Seguros No Vida/tablaOlhsson.xlsx")

Limpieza y estandarización de la base de datos

# Encuentra los NA en toda la base de datos
nas <- is.na(datos)

# Suma los NA por columna
nas_por_columna <- colSums(nas)

# Imprime la cantidad de NA por columna
print(nas_por_columna)
##        zona   clase veh    Edad veh       Bonus    duracion        nsin 
##           0           0           0           0           0           0 
##    costesin num.pólizas    freq.sin   severidad 
##           0           0           6         231
#Imputacion de datos faltantes con 0 para no tener NAs en la base de datos
#Donde se tiene NA se sustituye por 0 considerando que no se tuvo siniestros
datos$severidad <- ifelse(is.na(datos$severidad), 0, datos$severidad)
datos$freq.sin <- ifelse(is.na(datos$freq.sin), 0, datos$freq.sin)

Visualización inicial de las variables

Distribuciones de variables respuesta

# Distribución variable respuesta frecuencia
ggplot(data = datos, aes(x = freq.sin)) +
  geom_density(fill = "steelblue", alpha = 0.8) +
  geom_rug(alpha = 0.1) +
  scale_x_continuous(labels = scales::comma) +
  labs(title = "Distribución original") +
  theme_bw() 

# Distribución variable respuesta severidad
ggplot(data = datos, aes(x = severidad)) +
  geom_density(fill = "steelblue", alpha = 0.8) +
  geom_rug(alpha = 0.1) +
  scale_x_continuous(labels = scales::comma) +
  labs(title = "Distribución original") +
  theme_bw() 

Graficos de distribucion para cada variable numerica

variables_continuas <- datos %>% select(duracion, nsin, costesin, num.pólizas)

# Graficar la densidad de variables continuas
plot_list <- lapply(names(variables_continuas), function(variable_name) {
  ggplot(datos, aes(x = !!sym(variable_name))) + 
    geom_density() + 
    ggtitle(paste("Densidad de", variable_name))
})

# Mostrar los gráficos en una cuadrícula
gridExtra::grid.arrange(grobs = plot_list, ncol = 3)

#Matriz de correlacion entre las variables
# Instalar e cargar el paquete 'corrplot' si no lo has hecho aún
# install.packages("corrplot")
library(corrplot)
## corrplot 0.92 loaded
# Supongamos que 'df' es tu conjunto de datos
# Calcular la matriz de correlación
correlacion <- datos %>% select(duracion, nsin, costesin, num.pólizas, severidad, freq.sin)
correlation_matrix <- cor(correlacion)

# Imprimir la matriz de correlación
cat("\nMatriz de correlación:\n")
## 
## Matriz de correlación:
print(correlation_matrix)
##                duracion       nsin  costesin num.pólizas  severidad    freq.sin
## duracion     1.00000000 0.69470192 0.3736515  0.92701640 0.06514671 -0.05417345
## nsin         0.69470192 1.00000000 0.7163974  0.68881877 0.25588991  0.07441797
## costesin     0.37365152 0.71639739 1.0000000  0.33444135 0.62754705  0.11285339
## num.pólizas  0.92701640 0.68881877 0.3344413  1.00000000 0.07093177 -0.06329882
## severidad    0.06514671 0.25588991 0.6275470  0.07093177 1.00000000  0.17859102
## freq.sin    -0.05417345 0.07441797 0.1128534 -0.06329882 0.17859102  1.00000000
# Visualizar la matriz de correlación con un mapa de calor
# Puedes ajustar los parámetros según tus preferencias
corrplot(correlation_matrix, method = "color", 
         type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

# Titulo del gráfico
title("Matriz de correlación")

Existe poca correlación entre variables, por lo tanto no representa un riesgo para el funcionamiento de la red neuronal.

#Cambiar nombres
names(datos)[names(datos) == "clase veh"] <- "clase_vehiculo"
names(datos)[names(datos) == "Edad veh"] <- "Edad_vehiculo"

Preprocesado

Binarización de las variables categóricas para el correcto entrenamiento del modelo

Las variables zona, clase vehículo, edad vehículo y bonus, serán necesarias para la red neuronal y su predicción, sin embargo estas variables al ser categóricas, no pueden ser procesadas por la red. Para su correcta interpretación en el modelo, se crearán dummies con estas variables.

#Crear dummies para zona
zona_dummies <- model.matrix(~ factor(datos$zona) - 1)
colnames(zona_dummies) <- paste0("zona_", levels(factor(datos$zona)))
datos <- cbind(datos, zona_dummies)

#Crear dummies para clase_vehiculo
clase_vehiculo_dummies <- model.matrix(~ factor(datos$clase_vehiculo) - 1)
colnames(clase_vehiculo_dummies) <- paste0("clase_vehiculo_", levels(factor(datos$clase_vehiculo)))
datos <- cbind(datos, clase_vehiculo_dummies)

#Crear dummies para Edad_vehiculo
Edad_vehiculo_dummies <- model.matrix(~ factor(datos$Edad_vehiculo) - 1)
colnames(Edad_vehiculo_dummies) <- paste0("Edad_vehiculo_", levels(factor(datos$Edad_vehiculo)))
datos <- cbind(datos, Edad_vehiculo_dummies)

#Crear dummies para Bonus
Bonus_dummies <- model.matrix(~ factor(datos$Bonus) - 1)
colnames(Bonus_dummies) <- paste0("Bonus_", levels(factor(datos$Bonus)))
datos <- cbind(datos, Bonus_dummies)
#Eliminar las variables originales
datos <- datos[, !(names(datos) %in% c("zona", "clase_vehiculo", "Edad_vehiculo", "Bonus"))]

Estandarización de variables

Para las variables duración, nsin, costesin, num.polizas, freq.sin y severidad, se realizará una estandarización debido a que presentan diversas magnitudes, estas diferencias pueden influir ya que las variables que tengan una mayor escala o mayor varianza dominarán el modelo y podrían generar resultados sesgados.

\[ X_{\text{norm}} = \frac{X - \text{mean}(X)}{\text{sd}(X)} \]

#Debido a que las variables predictoras numéricas presentan diversas magnitudes, estas diferencias pueden influir ya que las variables que tengan una escala mayor o más varianza dominarán el modelo. 

# Calcula los parámetros de normalización
parametros_norm <- scale(datos)

Creación de modelos de redes

Separación de datos en conjunto de entrenamiento y prueba

# Dividir datos en conjuntos de entrenamiento y prueba
indice_entrenamiento <- sample(1:nrow(datos), 0.8 * nrow(datos))
conjunto_entrenamiento <- datos[indice_entrenamiento, ]
conjunto_prueba <- datos[-indice_entrenamiento, ]

#Modelos de redes neuronales

# Normaliza los datos de prueba solo para las variables seleccionadas
conjunto_entrenamiento_normalizado <- scale(conjunto_entrenamiento,
                                     center = colMeans(parametros_norm),
                                     scale = colSds(parametros_norm))
# Normaliza los datos de prueba solo para las variables seleccionadas
conjunto_prueba_normalizado <- scale(conjunto_prueba,
                                     center = colMeans(parametros_norm),
                                     scale = colSds(parametros_norm))
# Función para invertir la normalización
invertir_normalizacion <- function(predicciones_normalizadas, media, desviacion) {
  return (predicciones_normalizadas * desviacion + media)
}

Red neuronal para frecuencia

###Red neuronal de frecuencia -----------------------------------------------------
# Crear la red neuronal
modelo_frecuencia <- neuralnet(
    formula = freq.sin ~ .-freq.sin,
    data = conjunto_entrenamiento_normalizado,
    hidden = c(5, 3),
    linear.output = TRUE
  )
plot(modelo_frecuencia)

Una demostración gráfica del funcionamiento interno de la red neuronal donde se puede observar con las capas ocultas y nodos asignados.

# Realizar predicciones en el conjunto de prueba
predicciones_normalizadas_frecuencia <- predict(modelo_frecuencia, newdata = conjunto_prueba_normalizado)

# Invierte la normalización para obtener las predicciones originales
predicciones_originales_frecuencia <- invertir_normalizacion(predicciones_normalizadas_frecuencia,
                                                             mean(conjunto_prueba$freq.sin),
                                                             sd(conjunto_prueba$freq.sin))
# Calcular el error en el conjunto de prueba
error_frecuencia <- sqrt(mean((predicciones_originales_frecuencia - conjunto_prueba$freq.sin)^2))
cat("Error en el conjunto de prueba:", error_frecuencia, "\n")
## Error en el conjunto de prueba: 0.046644

Red neuronal para severidad

Para severidad se trabaja otro tipo de red neuronal utilizando cluster para el procesamiento de datos.

Esto permite reducir la carga en la computadora local, mejorando el tiempo de respuesta y al mismo tiempo aumentando la capacidad de cambios de parámetros para la red.

#Se realiza la binarización y estandarización de los datos para el entrenamiento y prueba del modelo.

transformer <- recipe(
                  formula = severidad ~ .,
                  data =  conjunto_entrenamiento
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

transformer
## 
## ── Recipe ──────────────────────────────────────────────────────────────────────
## 
## ── Inputs
## Number of variables by role
## outcome:    1
## predictor: 25
## 
## ── Operations
## • Removing rows with NA values in: all_predictors()
## • Sparse, unbalanced variable filter on: all_predictors()
## • Centering for: all_numeric(), -all_outcomes()
## • Scaling for: all_numeric(), -all_outcomes()
## • Dummy variables from: all_nominal(), -all_outcomes()
# Se entrena el objeto recipe
transformer_fit <- prep(transformer)

# Se aplican las transformaciones al conjunto de entrenamiento y de test
datos_train_prep <- bake(transformer_fit, new_data = conjunto_entrenamiento)
datos_test_prep  <- bake(transformer_fit, new_data = conjunto_prueba)
#Se inicia el cluster, el cual ayudará a reducir la carga en la computadora local y mejorará el tiempo de respuesta
h2o.init(nthreads = -1, max_mem_size = "4g")
## 
## H2O is not running yet, starting it now...
## 
## Note:  In case of errors look at the following log files:
##     C:\Users\Lenovo\AppData\Local\Temp\RtmpgJ4wQb\file5cf06d05505d/h2o_Lenovo_started_from_r.out
##     C:\Users\Lenovo\AppData\Local\Temp\RtmpgJ4wQb\file5cf03426178a/h2o_Lenovo_started_from_r.err
## 
## 
## Starting H2O JVM and connecting: . Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         4 seconds 529 milliseconds 
##     H2O cluster timezone:       Europe/Paris 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.44.0.3 
##     H2O cluster version age:    1 month and 11 days 
##     H2O cluster name:           H2O_started_from_R_Lenovo_uxr714 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   3.98 GB 
##     H2O cluster total cores:    8 
##     H2O cluster allowed cores:  8 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     R Version:                  R version 4.3.2 (2023-10-31 ucrt)
# Se eliminan los datos del cluster por si ya había sido iniciado.
h2o.removeAll()
h2o.no_progress()
#Se asignan los datos de entrenamiento y prueba al cluster
datos_train  <- as.h2o(conjunto_entrenamiento, key = "datos_train")
datos_test   <- as.h2o(conjunto_prueba, key = "datos_test")
# Espacio de búsqueda de cada hiperparámetro, donde se asignan las capas ocultas de la red neuronal
hiperparametros <- list(
                      epochs = c(50, 100, 500),
                      hidden = list(5, 10, 25, 50, c(10, 10))
                    )


# Búsqueda por validación cruzada
#Se asigna la variable respuesta y los predictores
variable_respuesta <- 'severidad'
predictores <- setdiff(colnames(datos_train), variable_respuesta)

#Se colocan los parámetros para la red neuronal
#Se incluye el algoritmo de deeplearning y validación cruzada

#Se está creando también un frame donde se evaluan varios modelos para lograr encontrar el que presente menor error.
grid <- h2o.grid(
              algorithm    = "deeplearning",
              activation   = "Rectifier",
              x            = predictores,
              y            = variable_respuesta,
              training_frame  = datos_train,
              nfolds       = 3, #validacion cruzada
              standardize  = FALSE,
              hyper_params = hiperparametros,
              search_criteria = list(strategy = "Cartesian"),
              seed         = 123,
              grid_id      = "grid"
          )
# Resultados de la red
resultados_grid <- h2o.getGrid(
                     sort_by = 'rmse',
                     grid_id = "grid",
                     decreasing = FALSE
                   )

data.frame(resultados_grid@summary_table)

En el grid de respuesta se puede evaluar la eficiencia de cada modelo evaluado por la validación cruzada del set de datos de entrenamiento y se toma el modelo evaluado con el menor error.

# Mejor modelo encontrado
modelo_final <- h2o.getModel(resultados_grid@model_ids[[1]])
#Se evaluan las predicciones para severidad con el set de datos de validación
predicciones <- h2o.predict(
                  object  = modelo_final,
                  newdata = datos_test
                )

predicciones <- predicciones %>%
                as_tibble() %>%
                mutate(valor_real = as.vector(datos_test$severidad))
# Calcular el error en el conjunto de prueba
cat("Error en el conjunto de prueba:", mean(resultados_grid@summary_table$rmse), "\n")
## Error en el conjunto de prueba: 19873.74

Conclusión

Resultados del error para ambos modelos de redes neuronales prediciendo los valores de número de frecuencias y severidad.

\[ RMSE = \sqrt{\frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{n}} \]

cat("RMSE para frecuencia:", error_frecuencia, "numero de siniestros", "\n")
## RMSE para frecuencia: 0.046644 numero de siniestros
cat("RMSE para severidad:", mean(resultados_grid@summary_table$rmse), "cantidad de euros", "\n")
## RMSE para severidad: 19873.74 cantidad de euros