1 ¿Por qué este material?

En las unidades anteriores trabajamos las metodologías (CRISP-DM, KDD, SEMMA, TDSP). En esta unidad vamos a recorrer las fases de CRISP-DM de punta a punta con un problema real del sector asegurador:

¿Cuánto le debe cobrar una aseguradora a cada asegurado de su cartera de autos?

Para responder vamos a usar una herramienta nueva: el Modelo Lineal Generalizado, que en todo el documento llamaremos GLM (por su sigla en inglés, Generalized Linear Model). No se preocupen por el nombre: es una “versión mejorada” de la regresión lineal que ya conocen, pensada para datos que no se comportan como una campana normal.

Vocabulario mínimo de seguros (lo usaremos todo el tiempo)

  • Póliza: el contrato de seguro de un cliente.
  • Siniestro: el evento que la aseguradora debe pagar (un choque, un robo).
  • Prima: lo que el cliente paga por el seguro.
  • Exposición: cuánto tiempo estuvo vigente la póliza. Una póliza de 6 meses tiene exposición 0,5 años.
  • Frecuencia: cuántos siniestros ocurren por cada año de exposición.
  • Severidad: cuánto cuesta, en promedio, cada siniestro.
  • Prima pura: lo que la aseguradora espera pagar en siniestros por un cliente en un año. Es frecuencia × severidad.

Cómo leer este documento. Las cajas grises que dicen “Para los curiosos” se pueden abrir con un clic. Tienen las fórmulas y los detalles técnicos. No son necesarias para entender el caso: si es su primer encuentro con estos temas, pueden saltarlas y volver después. Al final hay un glosario con todos los términos nuevos.

La ruta del documento es:

Parte Contenido Fase de la unidad
2 Del modelo lineal al GLM: la idea y las condiciones para usarlo Fundamentos
3 Dónde usar cada GLM: el mapa actuarial Fundamentos
4 Entendimiento del negocio y los datos III.1
5 Entender los datos III.2
6 Preparación de los datos III.3
7 Modelamiento: frecuencia, severidad y prima pura III.4
8 Evaluación, implementación y monitoreo III.5
9–10 Síntesis y glosario — 
# Si no los tienen instalados, quiten el # de las dos líneas siguientes y ejecútenlas una vez:
# install.packages(c("insuranceData", "dplyr", "ggplot2", "tidyr", "MASS", "statmod",
#                    "car", "kableExtra", "scales", "forcats", "gridExtra"))
library(insuranceData)  # la base de datos dataCar
library(dplyr)          # manipular datos
library(tidyr)          # reorganizar tablas
library(ggplot2)        # gráficos
library(MASS)           # modelo binomial negativo
library(statmod)        # residuales para diagnóstico y modelo Tweedie
library(car)            # revisar si las variables se repiten información (VIF)
library(kableExtra)     # tablas bonitas
library(scales)         # formatos de números y porcentajes

select <- dplyr::select  # MASS también tiene una función select(); usamos la de dplyr
theme_set(theme_minimal(base_size = 12))
tabla <- function(x, ...) {
  kbl(x, ...) |> kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
}

2 Del modelo lineal al GLM

2.1 Lo que ya sabemos: la regresión lineal

En la regresión lineal explicamos una variable \(Y\) con una recta (o un plano, si hay varias variables):

\[ Y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p + \text{error} \]

Detrás de esa ecuación hay tres condiciones que casi nunca decimos en voz alta:

  1. \(Y\) puede tomar cualquier valor, incluso negativos o con decimales.
  2. La dispersión es la misma para todos. Los datos se alejan de la recta más o menos lo mismo en cualquier punto.
  3. Los errores siguen una campana normal, simétrica alrededor de la recta.

Con datos de seguros, ninguna de las tres se cumple: el número de siniestros es un conteo (0, 1, 2…), nunca es negativo, la mayoría de clientes tiene cero y los costos tienen unos pocos valores enormes.

2.2 ¿Qué pasa si usamos una recta con conteos?

Simulemos el número de siniestros de 300 conductores según sus años de experiencia. Los conductores con menos experiencia tienen más siniestros.

set.seed(123)
n <- 300
experiencia <- runif(n, 0, 30)                                    # años de experiencia
siniestros  <- rpois(n, lambda = exp(1.2 - 0.12 * experiencia))   # conteos simulados
sim <- data.frame(experiencia, siniestros)

m_lm  <- lm(siniestros ~ experiencia, data = sim)                                # recta
m_glm <- glm(siniestros ~ experiencia, family = poisson(link = "log"), data = sim)  # GLM

grid <- data.frame(experiencia = seq(0, 40, by = 0.5))
grid$`Modelo lineal`  <- predict(m_lm, grid)
grid$`GLM Poisson` <- predict(m_glm, grid, type = "response")

grid |>
  pivot_longer(-experiencia, names_to = "modelo", values_to = "pred") |>
  ggplot(aes(experiencia, pred, color = modelo)) +
  geom_jitter(data = sim, aes(experiencia, siniestros), inherit.aes = FALSE,
              height = 0.1, alpha = 0.35) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Años de experiencia", y = "Número de siniestros",
       color = NULL, title = "La recta predice siniestros negativos") +
  theme(legend.position = "top")

range(predict(m_lm, grid))   # el valor mínimo que predice la recta
## [1] -1.468  2.356

La recta falla en tres cosas: predice un número de siniestros negativo para los conductores con mucha experiencia (¡imposible!), supone que todos se dispersan igual (pero los inexpertos varían mucho más) y obliga a una línea recta, cuando lo que se ve en los datos es una curva que baja rápido y luego se aplana. El GLM (la curva roja) resuelve las tres.

2.3 La idea del GLM en una frase

Un GLM es una regresión a la que le permitimos (1) escoger la distribución que se parece a nuestros datos y (2) trabajar en una escala distinta, para que sus predicciones siempre tengan sentido.

Para armar un GLM se toman tres decisiones:

Decisión 1. ¿Qué forma tienen los datos? (la distribución) Si \(Y\) es un conteo, usamos la distribución de Poisson; si es un monto de dinero siempre positivo y con valores muy grandes ocasionales, la Gamma; si es un sí/no, la Binomial. Es elegir el “molde” que mejor encaja con los datos.

Decisión 2. ¿Qué variables explican \(Y\)? (la parte lineal) Esto es igual que en la regresión de siempre: una suma de variables multiplicadas por coeficientes.

\[ \text{puntaje} = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \]

Decisión 3. ¿Cómo traducimos ese puntaje a la escala de \(Y\)? (la función de enlace) El puntaje puede ser cualquier número, incluso negativo. Pero un número de siniestros no puede ser negativo. El enlace es un traductor. El más usado en seguros es el logaritmo: el modelo trabaja con el puntaje y, para predecir, aplica la exponencial, que siempre da un número positivo.

puntajes <- c(-3, -1.5, 0, 1)
data.frame(puntaje = puntajes, `prediccion = exp(puntaje)` = exp(puntajes), check.names = FALSE) |>
  tabla(digits = 3, caption = "El enlace logarítmico convierte cualquier puntaje en un valor positivo")
El enlace logarítmico convierte cualquier puntaje en un valor positivo
puntaje prediccion = exp(puntaje)
-3.0 0.050
-1.5 0.223
0.0 1.000
1.0 2.718

Idea clave: la regresión lineal que ya conocen es un GLM: el que usa la distribución normal y no traduce nada (enlace “identidad”). Todo lo que aprendieron de regresión sigue sirviendo. El GLM solo agrega flexibilidad para los datos que no son “normales”.

Las distribuciones más usadas en un GLM
Distribución Tipo de Y Ejemplo en seguros ¿Cómo crece la variabilidad? Enlace usual
Normal Cualquier número Pocas veces aplica No crece: es constante Identidad
Binomial Sí / no (0 o 1) ¿Tuvo siniestro? Máxima cuando p = 0,5 Logit
Poisson Conteo: 0, 1, 2, … Número de siniestros Al ritmo de la media Logaritmo
Gamma Monto positivo y asimétrico Costo de un siniestro Más rápido que la media Logaritmo
Tweedie Muchos ceros + montos positivos Costo anual de una póliza Entre Poisson y Gamma Logaritmo
Para los curiosos: la familia exponencial y la función de varianza

Todas las distribuciones de la tabla pertenecen a la familia exponencial, cuya densidad se puede escribir como

\[ f(y;\theta,\phi) = \exp\left\{\frac{y\theta - b(\theta)}{\phi} + c(y,\phi)\right\}. \]

De esa forma se deduce que \(E(Y)=\mu=b'(\theta)\) y \(\operatorname{Var}(Y)=\phi\,V(\mu)\). La función de varianza \(V(\mu)\) es la “huella digital” de cada distribución: Normal \(V=1\); Poisson \(V=\mu\); Gamma \(V=\mu^2\); Gaussiana inversa \(V=\mu^3\); Tweedie \(V=\mu^p\) con \(1<p<2\); Binomial \(V=\mu(1-\mu)\).

Con notación formal, el GLM se escribe como \(g(\mu_i) = \eta_i = \beta_0 + \sum_j \beta_j x_{ij}\), donde \(\eta\) es el predictor lineal (el “puntaje”) y \(g\) es la función de enlace. El modelo fue propuesto por Nelder y Wedderburn (1972).

2.4 ¿Cómo encuentra el modelo los coeficientes?

En la regresión lineal usamos mínimos cuadrados: buscar la recta que deja los errores al cuadrado más pequeños. En un GLM la idea cambia un poco y se llama máxima verosimilitud. Suena complicado, pero la pregunta es sencilla:

¿Qué valores de los coeficientes hacen que los datos que observamos sean lo más “creíbles” posible?

Veámoslo con un ejemplo. En la simulación anterior, supongamos que todos los conductores tuvieran la misma frecuencia \(\lambda\). Probemos muchos valores de \(\lambda\) y calculemos qué tan probable sería observar nuestros datos con cada uno:

lambdas <- seq(0.3, 2.5, by = 0.01)
log_veros <- sapply(lambdas, function(l) sum(dpois(sim$siniestros, lambda = l, log = TRUE)))

mejor <- lambdas[which.max(log_veros)]

ggplot(data.frame(lambdas, log_veros), aes(lambdas, log_veros)) +
  geom_line(linewidth = 1, color = "#2c7fb8") +
  geom_vline(xintercept = mejor, linetype = "dashed", color = "#e6550d") +
  annotate("text", x = mejor + 0.05, y = min(log_veros) + 50, hjust = 0,
           label = paste0("Mejor valor: ", mejor), color = "#e6550d") +
  labs(x = "Valor propuesto para la frecuencia λ",
       y = "Qué tan creíbles son los datos\n(log-verosimilitud)",
       title = "Máxima verosimilitud: buscar la cima de la curva")

c(`λ que maximiza` = mejor, `Promedio de siniestros` = mean(sim$siniestros))
##         λ que maximiza Promedio de siniestros 
##                 0.9200                 0.9233

El valor que hace “más creíbles” los datos prácticamente coincide con el promedio de siniestros (la pequeña diferencia es porque probamos valores de 0,01 en 0,01), que es lo que la intuición diría. En un GLM, R hace exactamente esto, pero con todos los coeficientes a la vez. No hay que hacerlo a mano: la función glm() lo resuelve.

Para comparar modelos usaremos dos indicadores. Basta con saber cómo leerlos:

  • Devianza: mide qué tan lejos está el modelo de un ajuste perfecto. Es la “pariente” de la suma de errores al cuadrado de la regresión. Más pequeña, mejor.
  • AIC: es la devianza más un “castigo” por cada variable que agregamos, para no premiar modelos innecesariamente complicados. Más pequeño, mejor. Sirve para comparar modelos sobre los mismos datos.
Para los curiosos: devianza y algoritmo de estimación

La devianza se define como \(D = 2\,\phi\,[\ell(\text{modelo saturado}) - \ell(\text{modelo ajustado})]\), donde \(\ell\) es la log-verosimilitud y el modelo saturado es uno con un parámetro por observación (ajuste perfecto). La diferencia de devianzas entre dos modelos anidados da la prueba de razón de verosimilitud. El AIC es \(-2\ell + 2k\), con \(k\) el número de parámetros.

R maximiza la verosimilitud con el algoritmo IRLS (mínimos cuadrados ponderados iterativos): en cada paso resuelve una regresión ponderada y repite hasta que los coeficientes dejan de cambiar.

2.5 Cómo se interpretan los coeficientes: factores que multiplican

Con el enlace logarítmico, la predicción es

\[ \text{predicción} = e^{\beta_0}\times e^{\beta_1 x_1}\times\cdots\times e^{\beta_p x_p}. \]

Es decir, el efecto de cada variable multiplica en lugar de sumar. Por eso no leemos \(\beta\) directamente, sino \(e^{\beta}\), al que llamaremos factor:

  • Factor = 1,25 → esa característica aumenta la predicción un 25 %.
  • Factor = 0,80 → la disminuye un 20 %.
  • Factor = 1 → no cambia nada.

Ejemplo: si la frecuencia base es 0,15 siniestros al año, un conductor joven tiene factor 1,30 y su vehículo es costoso con factor 1,20, su frecuencia esperada es \(0{,}15 \times 1{,}30 \times 1{,}20 = 0{,}234\).

Así se construyen las tarifas de seguros en el mundo real: una tarifa base multiplicada por recargos y descuentos. Esta es la razón principal por la que los GLM son el estándar de la industria.

2.6 Las condiciones para usar un GLM (supuestos)

Antes de confiar en un modelo, hay que revisar que se cumplan ciertas condiciones. La tabla las presenta como preguntas que debemos hacernos.

Supuestos del GLM explicados como preguntas
Supuesto La pregunta que hay que hacerse Cómo lo revisamos
  1. Distribución adecuada
¿El molde (Poisson, Gamma…) se parece a la forma de mis datos? Gráficos de la variable: ¿cuántos ceros?, ¿cola larga?
  1. Variabilidad bien descrita
¿La variabilidad de los datos crece con la media como lo supone la distribución elegida? Gráfico de media contra varianza; indicador de dispersión cercano a 1
  1. Forma de la relación
Después de aplicar el logaritmo, ¿el efecto de cada variable es parecido a una recta? Gráfico de la frecuencia (en logaritmo) por grupos de la variable
  1. Independencia
¿Cada póliza es un caso distinto, o hay clientes repetidos que se influyen? Conocer cómo se construyó la base
  1. Exposición incluida
¿Tuve en cuenta que unas pólizas estuvieron vigentes más tiempo que otras? Incluir la exposición en el modelo (offset)
  1. Variables que no se repiten
¿Hay dos variables que dicen prácticamente lo mismo? Indicador VIF
  1. Sin datos que dominen el modelo
¿Hay unos pocos registros que, si los quito, cambian todo el modelo? Gráfico de residuales

Lo que cambia frente a la regresión lineal: ya no necesitamos que los datos sean normales ni que todos tengan la misma dispersión. Lo que se vuelve importante es elegir bien la distribución, para que el modelo describa correctamente cuánto varían los datos. Si nos equivocamos, el modelo puede acertar en los promedios, pero se equivoca en cuán seguro está de sus resultados, y entonces los p-valores dejan de ser confiables.

3 Dónde usar cada GLM: el mapa actuarial

La primera pregunta siempre es: ¿qué tipo de variable quiero predecir? La respuesta casi decide el modelo.

Qué GLM usar según la pregunta
Pregunta de negocio Tipo de Y Modelo Resultado
¿Esta póliza tendrá al menos un siniestro? Sí / no Binomial (logit) Probabilidad de siniestro
¿El cliente cancelará su póliza? Sí / no Binomial (logit) Tasa de cancelación
¿Qué probabilidad hay de que una persona fallezca este año? Sí / no Binomial (logit) Tabla de mortalidad
¿Cuántos siniestros tendrá un cliente por año? Conteo Poisson o Binomial negativa (log) Frecuencia
¿Cuántas muertes habrá en un grupo de edad? Conteo Poisson (log) Tasa de mortalidad
¿Cuánto cuesta un siniestro? Monto positivo con cola larga Gamma (log) Severidad (costo promedio)
¿Cuánto cuesta una póliza en el año (incluyendo los que no tienen siniestros)? Muchos ceros + montos Tweedie (log) Prima pura

En seguros de autos, la forma clásica de calcular lo que se debe cobrar separa el riesgo en dos piezas:

\[ \text{Prima pura} = \underbrace{\text{Frecuencia}}_{\text{¿cuántas veces?}} \times \underbrace{\text{Severidad}}_{\text{¿cuánto cuesta cada vez?}} \]

Así entendemos por qué un cliente es costoso: porque se accidenta más, porque sus accidentes son más caros, o por las dos cosas. Esta separación supone que el número de siniestros y su costo no dependen uno del otro, que es el supuesto habitual en los modelos actuariales.

4 Fase 1 · Entendimiento del negocio y los datos

4.1 El problema

Una aseguradora de autos quiere revisar su tarifa. Hoy cobra lo mismo a casi todos, y sospecha que:

  • algunos segmentos pagan de más y se van con la competencia;
  • otros pagan de menos y generan pérdidas.

Objetivo de negocio: cobrar a cada asegurado una prima acorde con su riesgo.

Objetivo analítico: estimar la prima pura anual (costo esperado de siniestros por año de exposición) para cada perfil de asegurado y, a partir de ella, la prima comercial.

Criterio de éxito: que, en datos que el modelo no ha visto, la prima pura total sea cercana a los siniestros reales (cociente real/esperado ≈ 1) y que el modelo ordene bien el riesgo: los asegurados que el modelo considera más caros deben ser, de verdad, los más costosos.

4.2 Los datos

Usamos dataCar del paquete insuranceData: pólizas de seguro de autos de un año (2004–2005), con 67 856 registros.

data(dataCar)
glimpse(dataCar)
## Rows: 67,856
## Columns: 11
## $ veh_value <dbl> 1.06, 1.03, 3.26, 4.14, 0.72, 2.01, 1.60, 1.47, 0.52, 0.38, …
## $ exposure  <dbl> 0.30390, 0.64887, 0.56947, 0.31759, 0.64887, 0.85421, 0.8542…
## $ clm       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, …
## $ numclaims <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, …
## $ claimcst0 <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, …
## $ veh_body  <fct> HBACK, HBACK, UTE, STNWG, HBACK, HDTOP, PANVN, HBACK, HBACK,…
## $ veh_age   <int> 3, 2, 2, 2, 4, 3, 3, 2, 4, 4, 2, 3, 2, 1, 3, 2, 3, 3, 4, 3, …
## $ gender    <fct> F, F, F, F, F, M, M, M, F, F, M, M, F, M, M, M, F, M, F, F, …
## $ area      <fct> C, A, E, D, C, C, A, B, A, B, A, C, C, A, B, C, F, C, D, C, …
## $ agecat    <int> 2, 4, 2, 2, 2, 4, 4, 6, 3, 4, 2, 4, 4, 5, 6, 4, 4, 4, 2, 3, …
## $ X_OBSTAT_ <fct> 01101    0    0    0, 01101    0    0    0, 01101    0    0 …
Diccionario de datos
Variable Descripción Rol en el modelo
veh_value Valor del vehículo (en decenas de miles) Predictor
exposure Fracción del año que la póliza estuvo vigente (0 a 1) Exposición (offset)
clm Indicador de al menos un siniestro (0/1) Respuesta binaria
numclaims Número de siniestros Respuesta de frecuencia
claimcst0 Costo total de los siniestros Respuesta de severidad
veh_body Tipo de carrocería Predictor
veh_age Edad del vehículo (1 = más nuevo, 4 = más viejo) Predictor
gender Sexo del conductor Predictor (¿ético?)
area Zona de residencia (A a F) Predictor
agecat Categoría de edad del conductor (1 = más joven, 6 = mayor) Predictor
X_OBSTAT_ Variable técnica sin información Se descarta

5 Fase 2 · Entender los datos

Esta fase decide la familia del GLM. Antes de modelar, miramos la forma de cada respuesta.

5.1 La frecuencia: muchos ceros y conteos pequeños

dataCar |>
  count(numclaims) |>
  mutate(porcentaje = percent(n / sum(n), accuracy = 0.01)) |>
  tabla(caption = "Distribución del número de siniestros")
Distribución del número de siniestros
numclaims n porcentaje
0 63232 93.19%
1 4333 6.39%
2 271 0.40%
3 18 0.03%
4 2 0.00%

Más del 93 % de las pólizas no tuvo siniestros. La respuesta es un conteo con media baja: candidata a Poisson.

5.2 La exposición: no todas las pólizas estuvieron vigentes el año completo

ggplot(dataCar, aes(exposure)) +
  geom_histogram(bins = 40, fill = "#2c7fb8", color = "white") +
  labs(x = "Exposición (fracción del año)", y = "Pólizas",
       title = "La exposición varía de unos días a un año completo")

Una póliza vigente 3 meses tiene menos oportunidad de tener un siniestro que una vigente 12 meses. Si no lo tenemos en cuenta, vamos a creer que los clientes de pólizas cortas son “más cuidadosos”, cuando en realidad solo estuvieron menos tiempo en la calle. Por eso no modelamos “siniestros por póliza”, sino siniestros por año de exposición. Más adelante veremos cómo se le dice esto al modelo con algo llamado offset.

5.3 La severidad: positiva y muy asimétrica

sev_obs <- dataCar |> filter(clm == 1) |> mutate(sev = claimcst0 / numclaims)
summary(sev_obs$sev)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     200     354     713    1916    1952   55922
ggplot(sev_obs, aes(sev)) +
  geom_histogram(bins = 50, fill = "#e6550d", color = "white") +
  scale_x_continuous(labels = comma) +
  labs(x = "Costo medio por siniestro", y = "Siniestros",
       title = "La severidad tiene una cola derecha muy larga")

La media es mucho mayor que la mediana: unos pocos siniestros muy costosos “jalan” el promedio hacia arriba. El histograma no se parece en nada a una campana normal; tiene una cola larga a la derecha. Para montos así, el molde candidato es la Gamma.

Pero hay algo más. Miremos los valores más frecuentes:

sev_obs |> count(claimcst0, sort = TRUE) |> head(5) |>
  mutate(porcentaje = percent(n / nrow(sev_obs), accuracy = 0.1)) |>
  tabla(caption = "Los cinco costos más repetidos")
Los cinco costos más repetidos
claimcst0 n porcentaje
200.0 695 15.0%
353.8 219 4.7%
389.9 94 2.0%
390.0 55 1.2%
367.7 42 0.9%

Un hallazgo para llevarle al negocio: alrededor del 15 % de los siniestros cuesta exactamente 200, y ninguno cuesta menos. Eso no es azar: puede ser un valor mínimo de registro, un costo fijo de atención o un efecto del deducible. Un analista no lo “limpia” por su cuenta; pregunta. Y, como veremos, este detalle afectará el ajuste del modelo de severidad.

5.4 Primeros patrones

frec_empirica <- function(var) {
  dataCar |>
    group_by(nivel = factor(.data[[var]])) |>
    summarise(polizas = n(), frecuencia = sum(numclaims) / sum(exposure), .groups = "drop") |>
    mutate(variable = var)
}

bind_rows(lapply(c("agecat", "veh_age", "area", "gender"), frec_empirica)) |>
  ggplot(aes(nivel, frecuencia)) +
  geom_col(fill = "#2c7fb8") +
  facet_wrap(~variable, scales = "free_x") +
  labs(x = NULL, y = "Siniestros por año de exposición",
       title = "Frecuencia observada por variable")

El patrón más claro: los conductores más jóvenes (agecat = 1) tienen mayor frecuencia.

6 Fase 3 · Preparación de los datos

6.1 Limpieza y transformación

autos <- dataCar |>
  filter(veh_value > 0) |>                       # vehículos sin valor registrado
  select(-X_OBSTAT_) |>
  mutate(
    agecat   = factor(agecat),
    veh_age  = factor(veh_age),
    area     = factor(area),
    gender   = factor(gender),
    # carrocerías con pocos datos se agrupan para estimar factores estables
    veh_body = forcats::fct_lump_min(veh_body, min = 1000, other_level = "OTRO"),
    veh_body = relevel(veh_body, ref = "SEDAN"),
    valor_cat = cut(veh_value, breaks = c(0, 1, 1.5, 2, 2.5, Inf), right = FALSE,
                    labels = c("<1", "1-1.5", "1.5-2", "2-2.5", "2.5+"))
  )

nrow(dataCar) - nrow(autos)   # registros eliminados
## [1] 53
table(autos$veh_body)
## 
## SEDAN HBACK HDTOP STNWG TRUCK   UTE  OTRO 
## 22232 18915  1579 16234  1742  4586  2515

¿Por qué convertir el valor del vehículo en tramos? El supuesto 3 pide que, después de aplicar el logaritmo, el efecto de la variable se parezca a una recta. El gráfico siguiente muestra que no es así: la frecuencia sube rápido para los vehículos baratos y luego se estabiliza. Una solución práctica, y muy común en seguros, es agrupar en tramos (menos de 1, de 1 a 1,5, etc.) y dejar que el modelo estime un factor para cada tramo.

autos |>
  mutate(tramo = ntile(veh_value, 20)) |>
  group_by(tramo) |>
  summarise(valor = median(veh_value), log_frec = log(sum(numclaims) / sum(exposure))) |>
  ggplot(aes(valor, log_frec)) +
  geom_point(size = 2.5, color = "#2c7fb8") +
  geom_smooth(method = "loess", se = FALSE, color = "#e6550d") +
  labs(x = "Valor del vehículo (mediana del tramo)", y = "log(frecuencia)",
       title = "Supuesto 3: la relación en escala log no es una recta")

6.2 Dos bases: frecuencia y severidad

set.seed(2026)
id_train <- sample(nrow(autos), size = 0.8 * nrow(autos))
train <- autos[id_train, ]
test  <- autos[-id_train, ]

# Severidad: solo pólizas con siniestro; costo medio por siniestro
train_sev <- train |> filter(clm == 1) |> mutate(sev = claimcst0 / numclaims)
test_sev  <- test  |> filter(clm == 1) |> mutate(sev = claimcst0 / numclaims)

data.frame(Base = c("Frecuencia (train)", "Frecuencia (test)", "Severidad (train)", "Severidad (test)"),
           Registros = c(nrow(train), nrow(test), nrow(train_sev), nrow(test_sev))) |>
  tabla()
Base Registros
Frecuencia (train) 54242
Frecuencia (test) 13561
Severidad (train) 3697
Severidad (test) 921

Separamos los datos en dos partes:

  • Entrenamiento (train, 80 %): con estos datos el modelo “estudia”.
  • Prueba (test, 20 %): se guardan y no se tocan hasta la evaluación. Funcionan como un examen con preguntas que el modelo nunca vio. Si le va bien ahí, podemos confiar en que funcionará con clientes nuevos.

Además, la severidad solo se puede estudiar en las pólizas que sí tuvieron siniestro, por eso creamos una base aparte para ella.

7 Fase 4 · Modelamiento

7.1 Verificando los supuestos antes de modelar

7.1.1 Relación media–varianza (supuesto 2)

Cada distribución tiene su propia “regla” sobre cómo crece la variabilidad cuando crece el promedio. En la Poisson, la varianza es igual a la media. En la Gamma, la varianza crece con el cuadrado de la media (mucho más rápido).

Para ver cuál regla siguen nuestros datos, formamos grupos de pólizas parecidas (misma edad, mismo tramo de valor, misma zona), calculamos la media y la varianza de cada grupo y las graficamos en escala logarítmica. En esa escala, la regla aparece como la pendiente de los puntos: pendiente 1 indica Poisson y pendiente 2 indica Gamma.

celdas_frec <- train |>
  group_by(agecat, valor_cat, area) |>
  summarise(media = mean(numclaims), varianza = var(numclaims), n = n(), .groups = "drop") |>
  filter(n >= 100)

celdas_sev <- train_sev |>
  group_by(agecat, area) |>
  summarise(media = mean(sev), varianza = var(sev), n = n(), .groups = "drop") |>
  filter(n >= 30)

p1 <- ggplot(celdas_frec, aes(log(media), log(varianza))) +
  geom_point(color = "#2c7fb8") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  labs(title = "Frecuencia", subtitle = "Línea: Var = μ (Poisson)")

ajuste_sev <- coef(lm(log(varianza) ~ log(media), data = celdas_sev))
p2 <- ggplot(celdas_sev, aes(log(media), log(varianza))) +
  geom_point(color = "#e6550d") +
  geom_smooth(method = "lm", se = FALSE, color = "grey40") +
  labs(title = "Severidad",
       subtitle = paste0("Pendiente estimada: ", round(ajuste_sev[2], 2), " (Gamma ⇒ 2)"))

gridExtra_ok <- requireNamespace("gridExtra", quietly = TRUE)
if (gridExtra_ok) gridExtra::grid.arrange(p1, p2, ncol = 2) else { print(p1); print(p2) }

  • Frecuencia: los puntos siguen de cerca la línea punteada (varianza = media). Poisson es razonable, aunque más adelante revisaremos si hay un poco más de variabilidad de la esperada.
  • Severidad: la varianza crece mucho más rápido que la media; la pendiente está por encima de 2. La Gamma es una buena primera opción, pero esta pendiente alta es una pista que retomaremos más adelante. Lo seguro es que la variabilidad no es constante, así que la regresión lineal clásica no serviría.

7.2 Modelo de frecuencia: Poisson con offset

Queremos predecir el número de siniestros con un modelo Poisson y enlace logarítmico. Pero hay que decirle al modelo cuánto tiempo estuvo vigente cada póliza. Eso se hace con un offset: una variable que entra al modelo sin coeficiente que estimar, solo para ajustar la escala. En la práctica, con el offset el modelo predice siniestros por año y luego los multiplica por el tiempo que la póliza estuvo vigente.

En R se escribe offset(log(exposure)). Usamos el logaritmo porque el modelo trabaja en escala logarítmica.

Para los curiosos: la ecuación con offset

\[ \log\left(E[N_i]\right) = \log(\text{exposición}_i) + \beta_0 + \sum_j \beta_j x_{ij} \;\;\Longrightarrow\;\; \frac{E[N_i]}{\text{exposición}_i} = e^{\beta_0 + \sum_j \beta_j x_{ij}} \]

El término \(\log(\text{exposición})\) entra con coeficiente fijo igual a 1.

Empezamos con todas las variables candidatas y usamos drop1(): esta función quita una variable a la vez y mide cuánto empeora el modelo. Si el p-valor (Pr(>Chi)) es pequeño (menor a 0,05), la variable sí aporta y conviene dejarla.

frec_completo <- glm(
  numclaims ~ agecat + area + veh_age + veh_body + valor_cat + gender + offset(log(exposure)),
  family = poisson(link = "log"), data = train
)
drop1(frec_completo, test = "Chisq")

La edad del conductor (agecat) aporta muchísimo; el tipo de carrocería (veh_body) y el valor del vehículo (valor_cat) también. La zona, la edad del vehículo y el sexo no muestran un aporte claro en la frecuencia. Nos quedamos con un modelo más sencillo, con las tres variables que sí aportan:

frec_pois <- glm(
  numclaims ~ agecat + veh_body + valor_cat + offset(log(exposure)),
  family = poisson(link = "log"), data = train
)
summary(frec_pois)
## 
## Call:
## glm(formula = numclaims ~ agecat + veh_body + valor_cat + offset(log(exposure)), 
##     family = poisson(link = "log"), data = train)
## 
## Coefficients:
##                Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)    -1.71930    0.06406  -26.84 < 0.0000000000000002 ***
## agecat2        -0.16701    0.06088   -2.74              0.00608 ** 
## agecat3        -0.21731    0.05928   -3.67              0.00025 ***
## agecat4        -0.25341    0.05924   -4.28    0.000018893182380 ***
## agecat5        -0.49752    0.06646   -7.49    0.000000000000071 ***
## agecat6        -0.49581    0.07651   -6.48    0.000000000091503 ***
## veh_bodyHBACK  -0.00511    0.04219   -0.12              0.90370    
## veh_bodyHDTOP  -0.10160    0.10685   -0.95              0.34165    
## veh_bodySTNWG  -0.08044    0.04809   -1.67              0.09438 .  
## veh_bodyTRUCK  -0.12793    0.10356   -1.24              0.21671    
## veh_bodyUTE    -0.22529    0.07316   -3.08              0.00207 ** 
## veh_bodyOTRO    0.16559    0.08052    2.06              0.03973 *  
## valor_cat1-1.5  0.11744    0.04727    2.48              0.01298 *  
## valor_cat1.5-2  0.14277    0.04936    2.89              0.00382 ** 
## valor_cat2-2.5  0.26156    0.06138    4.26    0.000020334181131 ***
## valor_cat2.5+   0.29570    0.05671    5.21    0.000000184340371 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 20315  on 54241  degrees of freedom
## Residual deviance: 20183  on 54226  degrees of freedom
## AIC: 27753
## 
## Number of Fisher Scoring iterations: 6

7.2.1 ¿Hay sobredispersión?

La Poisson tiene una regla estricta: la varianza debe ser igual a la media. Cuando los datos varían más de lo que permite esa regla, decimos que hay sobredispersión. Pasa mucho en seguros, porque hay clientes que, por razones que no vemos en los datos, son más propensos a accidentarse.

Para medirla usamos un indicador de dispersión: compara la variabilidad real con la que espera el modelo. Si el Poisson está bien, debe dar cerca de 1; si da bastante más de 1, hay sobredispersión.

dispersion_pearson <- sum(residuals(frec_pois, type = "pearson")^2) / df.residual(frec_pois)
dispersion_devianza <- deviance(frec_pois) / df.residual(frec_pois)
c(Pearson = dispersion_pearson, Devianza = dispersion_devianza)
##  Pearson Devianza 
##   1.3839   0.3722

Hay dos formas de calcular el indicador y aquí dan respuestas distintas. La versión basada en la devianza sale muy por debajo de 1 y parece decir que “todo está bien”, pero con tan pocos siniestros por póliza esa versión no es confiable. La versión de Pearson, que es la que hay que mirar en estos casos, está por encima de 1: hay sobredispersión moderada. Consecuencia práctica: el Poisson está demasiado seguro de sí mismo y algunas variables pueden parecer importantes sin serlo.

7.2.2 Alternativa: binomial negativa

La binomial negativa también sirve para conteos, pero es más flexible: permite que la varianza sea mayor que la media. Es como un Poisson con un “margen extra” de variabilidad, controlado por un parámetro llamado \(\theta\) (theta). En R se ajusta con glm.nb().

Para compararlos usamos el AIC (más pequeño es mejor) y miramos el error estándar (EE) de cada coeficiente, que indica qué tan seguro está el modelo de ese valor: un EE pequeño significa mucha seguridad.

Para los curiosos: la varianza de la binomial negativa \(\text{Var}(N)=\mu+\mu^2/\theta\). Cuando \(\theta\to\infty\) se recupera la Poisson. 
frec_nb <- glm.nb(
  numclaims ~ agecat + veh_body + valor_cat + offset(log(exposure)),
  data = train
)

data.frame(
  Modelo = c("Poisson", "Binomial negativa"),
  AIC = c(AIC(frec_pois), AIC(frec_nb)),
  theta = c(NA, frec_nb$theta)
) |> tabla(digits = 2, caption = "Comparación de modelos de frecuencia")
Comparación de modelos de frecuencia
Modelo AIC theta
Poisson 27753 NA
Binomial negativa 27731 2.61 
data.frame(
  Coeficiente = names(coef(frec_pois)),
  `Factor (Poisson)` = exp(coef(frec_pois)),
  `EE Poisson` = sqrt(diag(vcov(frec_pois))),
  `Factor (BN)` = exp(coef(frec_nb)),
  `EE BN` = sqrt(diag(vcov(frec_nb))),
  check.names = FALSE
) |> tabla(digits = 3, row.names = FALSE,
           caption = "Factores y errores estándar (EE) en ambos modelos")
Factores y errores estándar (EE) en ambos modelos
Coeficiente Factor (Poisson) EE Poisson Factor (BN) EE BN
(Intercept) 0.179 0.064 0.180 0.065
agecat2 0.846 0.061 0.844 0.062
agecat3 0.805 0.059 0.804 0.061
agecat4 0.776 0.059 0.775 0.060
agecat5 0.608 0.066 0.606 0.068
agecat6 0.609 0.077 0.607 0.078
veh_bodyHBACK 0.995 0.042 0.996 0.043
veh_bodyHDTOP 0.903 0.107 0.901 0.109
veh_bodySTNWG 0.923 0.048 0.923 0.049
veh_bodyTRUCK 0.880 0.104 0.879 0.106
veh_bodyUTE 0.798 0.073 0.798 0.074
veh_bodyOTRO 1.180 0.081 1.180 0.082
valor_cat1-1.5 1.125 0.047 1.122 0.048
valor_cat1.5-2 1.153 0.049 1.151 0.050
valor_cat2-2.5 1.299 0.061 1.295 0.063
valor_cat2.5+ 1.344 0.057 1.342 0.058

Lección: los factores de tarifa son prácticamente iguales en los dos modelos. Lo que cambia la sobredispersión es qué tan seguros podemos estar de esos factores. En esta cartera la diferencia en los errores estándar es pequeña, porque hay muy pocos siniestros por póliza y el “margen extra” casi no se nota. Aun así, la binomial negativa tiene menor AIC, así que la usaremos para la tarifa. En carteras con muchos siniestros por cliente (salud, flotas de empresas) la diferencia sí es grande.

7.2.3 ¿Y si olvidamos el offset?

frec_sin_offset <- glm.nb(numclaims ~ agecat + veh_body + valor_cat + veh_age, data = train)
frec_con_offset <- glm.nb(numclaims ~ agecat + veh_body + valor_cat + veh_age + offset(log(exposure)),
                          data = train)

data.frame(
  Término = c("Intercepto (frecuencia base)", "veh_age2", "veh_age3", "veh_age4"),
  `Con offset` = exp(coef(frec_con_offset))[c("(Intercept)", "veh_age2", "veh_age3", "veh_age4")],
  `Sin offset` = exp(coef(frec_sin_offset))[c("(Intercept)", "veh_age2", "veh_age3", "veh_age4")],
  check.names = FALSE
) |> tabla(digits = 3, row.names = FALSE)
Término Con offset Sin offset
Intercepto (frecuencia base) 0.180 0.069
veh_age2 1.108 1.223
veh_age3 0.990 1.126
veh_age4 1.003 1.185 
autos |> group_by(veh_age) |> summarise(exposicion_media = mean(exposure)) |> tabla(digits = 3)
veh_age exposicion_media
1 0.436
2 0.478
3 0.476
4 0.474

Olvidar el offset produce dos errores.

  1. La frecuencia base (el intercepto) queda en “siniestros por póliza” y no “por año”, y sale en menos de la mitad. Si cobramos con ese número, la aseguradora cobraría mucho menos de lo necesario y perdería dinero.
  2. Los vehículos nuevos (veh_age = 1) estuvieron asegurados menos tiempo en promedio. Sin offset parecen más seguros de lo que son, y a los vehículos viejos se les cobraría un recargo injusto. Con el offset, ese supuesto efecto de la edad del vehículo prácticamente desaparece.

7.2.4 Diagnóstico con residuales cuantílicos

En regresión lineal revisábamos los residuales (lo que el modelo no logra explicar) para ver si el modelo estaba bien. Con conteos, los residuales comunes forman “rayas” que no se pueden leer. Por eso usamos los residuales cuantílicos: una transformación que convierte los residuales en valores que, si el modelo es correcto, se comportan como una normal.

Se leen con un gráfico QQ, como en regresión lineal: si los puntos caen sobre la recta, el modelo describe bien los datos; si se alejan, algo falla.

Para los curiosos: residuales cuantílicos aleatorizados Propuestos por Dunn y Smyth (1996). Para cada observación se toma la probabilidad acumulada \(F(y_i;\hat\mu_i)\) del modelo ajustado (con un valor aleatorio dentro del salto si la variable es discreta) y se transforma con la inversa de la normal estándar: \(r_i=\Phi^{-1}(u_i)\). Si el modelo es correcto, \(r_i\sim N(0,1)\). 
set.seed(1)
rq_pois <- qresiduals(frec_pois)
rq_nb   <- qresiduals(frec_nb)

data.frame(Poisson = rq_pois, `Binomial negativa` = rq_nb, check.names = FALSE) |>
  pivot_longer(everything(), names_to = "modelo", values_to = "residual") |>
  ggplot(aes(sample = residual)) +
  stat_qq(alpha = 0.3, size = 0.6, color = "#2c7fb8") +
  stat_qq_line(color = "#e6550d") +
  facet_wrap(~modelo) +
  labs(title = "QQ-plot de residuales cuantílicos aleatorizados",
       x = "Cuantiles teóricos N(0,1)", y = "Residual cuantílico")

En los dos modelos los puntos siguen la recta casi a la perfección: el modelo de frecuencia describe bien los datos.

7.2.5 ¿Hay variables que se repiten? (multicolinealidad)

Si dos variables dicen prácticamente lo mismo (por ejemplo, el valor del vehículo y su precio de compra), el modelo no sabe a cuál atribuirle el efecto y los factores se vuelven inestables. El indicador VIF lo detecta.

vif(frec_nb)
##            GVIF Df GVIF^(1/(2*Df))
## agecat    1.042  5           1.004
## veh_body  1.518  6           1.035
## valor_cat 1.474  4           1.050

Miramos la última columna: valores cercanos a 1 indican que no hay problema; valores mayores a 2 o 3 serían una señal de alerta. Aquí todo está bien.

7.2.6 Interpretación: factores de tarifa

Convertimos cada coeficiente en su factor (\(e^{\beta}\)) y calculamos un intervalo de confianza del 95 %: el rango donde razonablemente está el verdadero valor del factor. Si el intervalo no incluye el 1, la característica tiene un efecto claro.

factores_frec <- exp(cbind(Factor = coef(frec_nb), confint.default(frec_nb)))
factores_frec |> tabla(digits = 3, caption = "Factores multiplicativos de frecuencia (binomial negativa)")
Factores multiplicativos de frecuencia (binomial negativa)
Factor 2.5 % 97.5 %
(Intercept) 0.180 0.158 0.205
agecat2 0.844 0.747 0.953
agecat3 0.804 0.714 0.905
agecat4 0.775 0.688 0.872
agecat5 0.606 0.531 0.692
agecat6 0.607 0.521 0.707
veh_bodyHBACK 0.996 0.916 1.084
veh_bodyHDTOP 0.901 0.728 1.116
veh_bodySTNWG 0.923 0.838 1.016
veh_bodyTRUCK 0.879 0.714 1.081
veh_bodyUTE 0.798 0.689 0.923
veh_bodyOTRO 1.180 1.004 1.387
valor_cat1-1.5 1.122 1.021 1.233
valor_cat1.5-2 1.151 1.043 1.270
valor_cat2-2.5 1.295 1.145 1.464
valor_cat2.5+ 1.342 1.198 1.503 
data.frame(termino = rownames(factores_frec), factores_frec, check.names = FALSE) |>
  filter(termino != "(Intercept)") |>
  ggplot(aes(x = Factor, y = reorder(termino, Factor))) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  geom_errorbarh(aes(xmin = `2.5 %`, xmax = `97.5 %`), height = 0.2, color = "grey50") +
  geom_point(size = 3, color = "#2c7fb8") +
  labs(x = "Factor multiplicativo sobre la frecuencia", y = NULL,
       title = "Factores con intervalo de confianza del 95 %",
       subtitle = "A la derecha de 1: recargo · a la izquierda de 1: descuento")

Lectura para la gerencia, por ejemplo: un conductor de la categoría de edad 6 tiene una frecuencia de siniestros 39 % menor que un conductor de la categoría 1 (la más joven), con las demás características iguales.

7.3 Modelo de severidad: de la Gamma a la Gaussiana inversa

Ahora modelamos cuánto cuesta cada siniestro, usando solo las pólizas que tuvieron al menos uno. Le damos más peso (weights) a las pólizas con varios siniestros, porque una póliza con 3 siniestros nos dice más sobre el costo típico que una con 1.

De nuevo usamos drop1() para ver qué variables aportan. Aquí la prueba se llama F en lugar de Chi-cuadrado, porque la Gamma tiene un parámetro extra de dispersión que también se estima; la lectura es la misma: p-valor pequeño = la variable aporta.

sev_completo <- glm(
  sev ~ agecat + area + veh_age + veh_body + valor_cat + gender,
  family = Gamma(link = "log"), weights = numclaims, data = train_sev
)
drop1(sev_completo, test = "F")

En la severidad aportan con claridad la edad del conductor, la zona y el sexo. El tipo de carrocería tiene un aporte débil y, para mantener el modelo sencillo, lo dejamos por fuera. Curiosamente, son variables distintas de las que importaban en la frecuencia: lo que hace que alguien se accidente más no es lo mismo que hace que sus accidentes cuesten más.

sev_gamma <- glm(
  sev ~ agecat + area + gender,
  family = Gamma(link = "log"), weights = numclaims, data = train_sev
)
summary(sev_gamma)
## 
## Call:
## glm(formula = sev ~ agecat + area + gender, family = Gamma(link = "log"), 
##     data = train_sev, weights = numclaims)
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   7.6133     0.1044   72.90 <0.0000000000000002 ***
## agecat2      -0.1646     0.1083   -1.52              0.1288    
## agecat3      -0.2504     0.1051   -2.38              0.0173 *  
## agecat4      -0.2548     0.1053   -2.42              0.0156 *  
## agecat5      -0.3503     0.1181   -2.97              0.0030 ** 
## agecat6      -0.2358     0.1358   -1.74              0.0826 .  
## areaB        -0.0288     0.0852   -0.34              0.7354    
## areaC         0.0987     0.0779    1.27              0.2052    
## areaD         0.0400     0.1054    0.38              0.7042    
## areaE         0.1352     0.1128    1.20              0.2310    
## areaF         0.3792     0.1305    2.90              0.0037 ** 
## genderM       0.1741     0.0577    3.02              0.0026 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 3.183)
## 
##     Null deviance: 6074.8  on 3696  degrees of freedom
## Residual deviance: 5967.2  on 3685  degrees of freedom
## AIC: 67118
## 
## Number of Fisher Scoring iterations: 7

El parámetro de dispersión (Dispersion parameter en la salida) indica qué tanto varían los costos. Su raíz cuadrada es el coeficiente de variación: 1.78. Eso significa que la desviación estándar del costo es casi el doble de su promedio: los costos de los siniestros son muy variables, como ya anticipaba el histograma.

7.3.1 ¿Se cumple el supuesto de la distribución?

Dos pistas ya nos advertían: la pendiente media–varianza de la severidad salió por encima de 2, y el 15 % de los costos está pegado en 200. Revisemos los residuales cuantílicos.

set.seed(1)
rq_gamma <- qresiduals(sev_gamma)
ggplot(data.frame(rq = rq_gamma), aes(sample = rq)) +
  stat_qq(alpha = 0.4, color = "#e6550d") +
  stat_qq_line() +
  labs(title = "Severidad Gamma: residuales cuantílicos",
       x = "Cuantiles teóricos N(0,1)", y = "Residual cuantílico")

La Gamma no ajusta bien: los puntos se alejan de la recta en los dos extremos.

  • A la izquierda, los puntos quedan aplastados en una línea horizontal: son los siniestros de exactamente 200, que forman un “piso”.
  • A la derecha, los puntos se disparan por encima de la recta: hay siniestros más costosos de los que la Gamma considera posibles.

Esto es justo lo que pasa en la vida real: el primer modelo casi nunca es el último. El diagnóstico nos dice que debemos probar otra cosa.

7.3.2 Alternativa: Gaussiana inversa

Recordemos la pista: la pendiente media–varianza de la severidad era mayor que 2, es decir, la variabilidad crece todavía más rápido de lo que supone la Gamma. Existe otra distribución para montos positivos que admite eso: la Gaussiana inversa. Tiene una cola derecha más larga, así que acepta mejor los siniestros muy costosos. En R basta con cambiar family = Gamma(...) por family = inverse.gaussian(...).

sev_ig <- glm(
  sev ~ agecat + area + gender,
  family = inverse.gaussian(link = "log"), weights = numclaims, data = train_sev,
  control = glm.control(maxit = 100)
)

set.seed(1)
rq_ig <- qresiduals(sev_ig)

data.frame(Gamma = rq_gamma, `Gaussiana inversa` = rq_ig, check.names = FALSE) |>
  pivot_longer(everything(), names_to = "modelo", values_to = "residual") |>
  ggplot(aes(sample = residual)) +
  stat_qq(alpha = 0.4, color = "#e6550d") +
  stat_qq_line() +
  facet_wrap(~modelo) +
  labs(title = "Residuales cuantílicos: Gamma vs. Gaussiana inversa",
       x = "Cuantiles teóricos N(0,1)", y = "Residual cuantílico")

data.frame(Modelo = c("Gamma", "Gaussiana inversa"),
           `V(μ)` = c("μ²", "μ³"),
           AIC = c(AIC(sev_gamma), AIC(sev_ig)), check.names = FALSE) |>
  tabla(digits = 0, caption = "Comparación de modelos de severidad")
Comparación de modelos de severidad
Modelo V(μ) AIC
Gamma μ² 67118
Gaussiana inversa μ³ 65224

Decisión: la Gaussiana inversa tiene un AIC mucho menor y sus puntos siguen la recta en el centro y a la derecha. Para la tarifa usaremos la Gaussiana inversa.

El “piso” en 200 sigue ahí: ninguna distribución lo resuelve, porque es una característica de cómo se registran los datos. En la práctica se trataría aparte, por ejemplo modelando solo lo que cuesta cada siniestro por encima de 200. Algo parecido se hace con los siniestros gigantes: se analizan por separado y su costo se reparte entre todos como un pequeño recargo.

sev_final <- sev_ig

exp(cbind(Factor = coef(sev_final), confint.default(sev_final))) |>
  tabla(digits = 3, caption = "Factores multiplicativos de severidad (Gaussiana inversa)")
Factores multiplicativos de severidad (Gaussiana inversa)
Factor 2.5 % 97.5 %
(Intercept) 2041.347 1638.301 2543.547
agecat2 0.848 0.673 1.068
agecat3 0.783 0.626 0.978
agecat4 0.772 0.618 0.964
agecat5 0.704 0.553 0.896
agecat6 0.772 0.585 1.018
areaB 0.966 0.822 1.135
areaC 1.105 0.950 1.286
areaD 1.033 0.844 1.265
areaE 1.140 0.910 1.429
areaF 1.450 1.083 1.941
genderM 1.179 1.051 1.322

Una pregunta que no es estadística: la variable gender resulta significativa en la severidad. ¿Es aceptable cobrar distinto por sexo? En la Unión Europea, por ejemplo, la tarificación por sexo en seguros está prohibida desde 2012. Que una variable sea significativa no la vuelve legítima: la fase de entendimiento del negocio también incluye el marco regulatorio y ético.

7.4 La prima pura

Ya tenemos las dos piezas. Para cada cliente calculamos:

\[ \text{Prima pura} = \text{Siniestros esperados en un año} \times \text{Costo esperado de cada siniestro} \]

Para obtener los siniestros en un año completo, le pedimos al modelo de frecuencia la predicción con exposición igual a 1.

prima_pura <- function(datos) {
  anual <- mutate(datos, exposure = 1)
  frecuencia <- predict(frec_nb, newdata = anual, type = "response")
  severidad  <- predict(sev_final, newdata = anual, type = "response")
  data.frame(frecuencia, severidad, prima_pura = frecuencia * severidad)
}

train <- bind_cols(train, prima_pura(train))
test  <- bind_cols(test,  prima_pura(test))

summary(test$prima_pura)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     121     234     281     296     337     843
ggplot(test, aes(prima_pura)) +
  geom_histogram(bins = 50, fill = "#31a354", color = "white") +
  labs(x = "Prima pura anual estimada", y = "Pólizas",
       title = "Distribución de la prima pura en la base de prueba")

7.4.1 Una alternativa directa: Tweedie

¿Y si en lugar de dos modelos usamos uno solo? La distribución Tweedie está diseñada para variables que tienen muchos ceros y, cuando no son cero, montos positivos: exactamente como el costo anual de una póliza. Con ella se modela la prima pura directamente.

Para los curiosos: la Tweedie Con \(1<p<2\), la Tweedie es una distribución Poisson compuesta con sumandos Gamma: un número Poisson de siniestros, cada uno con costo Gamma. Su función de varianza es \(V(\mu)=\mu^p\). En statmod::tweedie(), var.power es \(p\) y link.power = 0 indica enlace logarítmico. 
tweedie_pp <- glm(
  I(claimcst0 / exposure) ~ agecat + veh_body + valor_cat + area + gender,
  family = tweedie(var.power = 1.6, link.power = 0),   # link.power = 0 ⇒ enlace log
  weights = exposure, data = train
)
test$prima_tweedie <- predict(tweedie_pp, newdata = test, type = "response")

cor(test$prima_pura, test$prima_tweedie)
## [1] 0.9174

Ambos enfoques producen primas similares: la correlación entre ellas es alta (1 sería idénticas). La ventaja del enfoque frecuencia × severidad es que explica de dónde viene el riesgo; la del Tweedie, que es un solo modelo.

8 Fase 5 · Evaluación, implementación y monitoreo

8.1 Evaluación en datos no vistos

Llegó el momento del “examen”: usamos la base de prueba, que el modelo nunca vio.

8.1.1 ¿El total cuadra?

Sumamos lo que el modelo esperaba pagar en siniestros y lo comparamos con lo que realmente se pagó. El indicador real / esperado se lee así: 1 significa que acertó; 1,10 que los siniestros fueron un 10 % más de lo previsto (la aseguradora cobró de menos); 0,90 que fueron un 10 % menos (cobró de más).

test <- test |> mutate(esperado = prima_pura * exposure)  # costo esperado según exposición real

eval_total <- test |>
  summarise(`Siniestros reales` = sum(claimcst0),
            `Siniestros esperados` = sum(esperado)) |>
  mutate(`Real / esperado` = `Siniestros reales` / `Siniestros esperados`)

eval_total |> tabla(digits = 3, format.args = list(big.mark = ","))
Siniestros reales Siniestros esperados Real / esperado
1,859,359 1,851,337 1.004

8.1.2 ¿Cuadra por segmento?

test |>
  group_by(agecat) |>
  summarise(polizas = n(), siniestros = sum(numclaims),
            real = sum(claimcst0), esperado = sum(esperado)) |>
  mutate(`Real / esperado` = real / esperado) |>
  tabla(digits = 3, format.args = list(big.mark = ","),
        caption = "Cociente real/esperado por categoría de edad")
Cociente real/esperado por categoría de edad
agecat polizas siniestros real esperado Real / esperado
1 1,194 112 336,168 259,648 1.295
2 2,602 200 393,350 401,065 0.981
3 3,154 223 387,067 437,934 0.884
4 3,257 234 425,449 425,966 0.999
5 2,060 133 196,946 192,332 1.024
6 1,294 89 120,380 134,392 0.896

Un cociente cercano a 1 en cada segmento indica que el modelo no está subsidiando a unos grupos con otros. Ojo con la lectura: cada segmento tiene solo unas decenas o pocos cientos de siniestros y la severidad es muy variable, así que desviaciones de 10 % a 30 % pueden ser ruido. Antes de cambiar la tarifa de un segmento, se confirma con más periodos o con intervalos de confianza.

8.1.3 ¿Ordena bien el riesgo? Gráfico de lift

Que el total cuadre no basta: un modelo que le cobre lo mismo a todos también puede cuadrar en total. Lo valioso es que distinga a los clientes baratos de los caros. Para comprobarlo, ordenamos las pólizas de prueba de menor a mayor prima estimada, las partimos en 10 grupos del mismo tamaño (deciles) y comparamos en cada grupo el costo que el modelo predijo con el que realmente ocurrió. A este gráfico se le llama lift.

lift <- test |>
  mutate(decil = ntile(prima_pura, 10)) |>
  group_by(decil) |>
  summarise(Predicho = sum(esperado) / sum(exposure),
            Observado = sum(claimcst0) / sum(exposure), .groups = "drop")

lift |>
  pivot_longer(-decil, names_to = "tipo", values_to = "costo") |>
  ggplot(aes(factor(decil), costo, color = tipo, group = tipo)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  scale_color_manual(values = c(Observado = "#e6550d", Predicho = "#2c7fb8")) +
  labs(x = "Decil de prima pura estimada (1 = menor riesgo)",
       y = "Costo por año de exposición", color = NULL,
       title = "Lift: ¿el modelo separa a los asegurados baratos de los caros?") +
  theme(legend.position = "top")

with(lift, Observado[10] / Observado[1])
## [1] 2.77

Si el modelo no sirviera para distinguir clientes, la línea naranja (lo observado) sería plana. En cambio, sube: los clientes que el modelo considera más riesgosos de verdad cuestan más. El número de abajo indica cuántas veces más cuesta el grupo más riesgoso frente al menos riesgoso.

Nota honesta: con variables solo del vehículo y de la edad, el poder de discriminación es limitado. En la industria se agregan historial de siniestros, kilometraje, uso del vehículo o datos telemáticos. El modelo es tan bueno como los datos que lo alimentan.

8.2 Implementación: de la prima pura a la prima comercial

La prima pura alcanza solo para pagar los siniestros esperados. Pero la aseguradora también tiene gastos (salarios, comisiones a los asesores, tecnología) y necesita una ganancia. La prima comercial, que es el precio que paga el cliente, los incluye:

\[ \text{Prima comercial} = \frac{\text{Prima pura}}{1 - g - u} \]

donde \(g\) es la proporción de gastos (administración, comisiones) y \(u\) el margen de utilidad sobre la prima. Por ejemplo, con gastos del 25 % y utilidad del 5 %, una prima pura de 300 se convierte en una prima comercial de \(300 / 0{,}70 \approx 429\).

cotizar <- function(agecat, veh_body, valor, area, gender,
                    gastos = 0.25, utilidad = 0.05) {
  perfil <- data.frame(
    agecat   = factor(agecat, levels = levels(autos$agecat)),
    veh_body = factor(veh_body, levels = levels(autos$veh_body)),
    valor_cat = cut(valor, breaks = c(0, 1, 1.5, 2, 2.5, Inf), right = FALSE,
                    labels = levels(autos$valor_cat)),
    area     = factor(area, levels = levels(autos$area)),
    gender   = factor(gender, levels = levels(autos$gender)),
    exposure = 1
  )
  pp <- prima_pura(perfil)
  data.frame(pp, prima_comercial = pp$prima_pura / (1 - gastos - utilidad))
}

# Dos perfiles para comparar
bind_rows(
  `Joven, hatchback, zona C` = cotizar(1, "HBACK", 1.2, "C", "M"),
  `Mayor, sedán, zona A`     = cotizar(6, "SEDAN", 1.2, "A", "F"),
  .id = "Perfil"
) |> tabla(digits = 3, row.names = FALSE, caption = "Cotización para dos perfiles")
Cotización para dos perfiles
Perfil frecuencia severidad prima_pura prima_comercial
Joven, hatchback, zona C 0.201 2659 535.4 764.9
Mayor, sedán, zona A 0.123 1576 193.4 276.3

En producción, esta función se convertiría en un servicio (una API, un aplicativo Shiny o una tabla de factores cargada en el sistema de emisión de pólizas).

8.2.1 La tarifa como tabla de factores

base <- exp(coef(frec_nb)["(Intercept)"]) * exp(coef(sev_final)["(Intercept)"])
cat("Prima pura base (perfil de referencia):", round(base, 1), "\n")
## Prima pura base (perfil de referencia): 367.7
rel_frec <- exp(coef(frec_nb)[-1])
rel_sev  <- exp(coef(sev_final)[-1])
data.frame(Término = c(names(rel_frec), names(rel_sev)),
           Componente = c(rep("Frecuencia", length(rel_frec)), rep("Severidad", length(rel_sev))),
           Factor = c(rel_frec, rel_sev)) |>
  tabla(digits = 3, row.names = FALSE, caption = "Tarifa expresada como factores multiplicativos")
Tarifa expresada como factores multiplicativos
Término Componente Factor
agecat2 Frecuencia 0.844
agecat3 Frecuencia 0.804
agecat4 Frecuencia 0.775
agecat5 Frecuencia 0.606
agecat6 Frecuencia 0.607
veh_bodyHBACK Frecuencia 0.996
veh_bodyHDTOP Frecuencia 0.901
veh_bodySTNWG Frecuencia 0.923
veh_bodyTRUCK Frecuencia 0.879
veh_bodyUTE Frecuencia 0.798
veh_bodyOTRO Frecuencia 1.180
valor_cat1-1.5 Frecuencia 1.122
valor_cat1.5-2 Frecuencia 1.151
valor_cat2-2.5 Frecuencia 1.295
valor_cat2.5+ Frecuencia 1.342
agecat2 Severidad 0.848
agecat3 Severidad 0.783
agecat4 Severidad 0.772
agecat5 Severidad 0.704
agecat6 Severidad 0.772
areaB Severidad 0.966
areaC Severidad 1.105
areaD Severidad 1.033
areaE Severidad 1.140
areaF Severidad 1.450
genderM Severidad 1.179

La prima de cualquier asegurado se obtiene multiplicando la prima base por los factores que le corresponden. Así se comunica una tarifa a un área comercial o a un regulador.

8.3 Monitoreo: el modelo envejece

Un modelo implementado no se deja solo. Los costos de reparación suben, cambia la mezcla de clientes y aparecen nuevos riesgos. Simulemos 12 meses del año siguiente en los que la inflación de los costos de reparación es del 1,5 % mensual, mientras la tarifa sigue igual.

Como la severidad es muy volátil, un mes aislado dice poco. Por eso se monitorea el cociente acumulado en lo corrido del año: cada mes suma más información y la banda de variación normal se estrecha.

set.seed(2026)
polizas_mes <- 4000

# Una cartera mensual simulada: pólizas remuestreadas del test, con inflación
simular_mes <- function(inflacion) {
  muestra <- test[sample(nrow(test), polizas_mes, replace = TRUE), ]
  c(real = sum(muestra$claimcst0) * inflacion, esperado = sum(muestra$esperado))
}

# Escenario sin cambios: 300 años simulados para construir las bandas de control
sin_cambio <- replicate(300, {
  meses <- sapply(1:12, function(m) simular_mes(1))
  cumsum(meses["real", ]) / cumsum(meses["esperado", ])
})
bandas <- data.frame(mes = 1:12,
                     inf = apply(sin_cambio, 1, quantile, 0.05),
                     sup = apply(sin_cambio, 1, quantile, 0.95))

# Año observado: con inflación de 1,5 % mensual
meses_obs <- sapply(1:12, function(m) simular_mes(1.015^m))
monitoreo <- bandas |>
  mutate(acumulado = cumsum(meses_obs["real", ]) / cumsum(meses_obs["esperado", ]),
         estado = ifelse(acumulado > sup, "Alerta", "Dentro de la banda"))

ggplot(monitoreo, aes(mes)) +
  geom_ribbon(aes(ymin = inf, ymax = sup), fill = "#31a354", alpha = 0.15) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_line(aes(y = acumulado), color = "grey40") +
  geom_point(aes(y = acumulado, color = estado), size = 3) +
  scale_color_manual(values = c(`Dentro de la banda` = "#31a354", Alerta = "#e6550d"), name = NULL) +
  scale_x_continuous(breaks = 1:12) +
  labs(x = "Mes del año siguiente", y = "Real / esperado (acumulado)",
       title = "Tablero de monitoreo del modelo",
       subtitle = "Banda verde: variación normal sin cambios (percentiles 5–95)") +
  theme(legend.position = "top")

Al principio del año la banda es ancha y la inflación se confunde con el ruido. A medida que se acumulan meses, la banda se estrecha y la desviación sale a la luz.

¿Qué hacer cuando se dispara la alerta? Primero se diagnostica: ¿es un mes atípico o una tendencia? Si es tendencia, se distingue entre un cambio de nivel (inflación: basta con ajustar la prima base) y un cambio de estructura (cambió el riesgo relativo entre segmentos: hay que reentrenar el modelo). Así el ciclo CRISP-DM vuelve a empezar.

9 Síntesis

Fase Lo que hicimos Supuesto o concepto clave
Negocio Definir la prima pura como objetivo y el cociente real/esperado como criterio de éxito Independencia frecuencia–severidad; marco ético
Entender los datos Identificar ceros, exposición variable y asimetría: eso decidió las familias Familia exponencial (supuesto 1)
Preparación Agrupar categorías, tramos de valor, separar frecuencia y severidad, train/test Linealidad en la escala del enlace (supuesto 3)
Modelamiento Poisson → sobredispersión → Binomial negativa; Gamma → falla el diagnóstico → Gaussiana inversa para severidad; Tweedie como alternativa Relación media–varianza (supuesto 2), offset (supuesto 5)
Evaluación Real/esperado total y por segmento; lift por deciles Validación fuera de muestra
Implementación Función cotizar() y tarifa como tabla de factores Prima comercial = PP / (1 − g − u)
Monitoreo Tablero real/esperado con bandas de control Deriva de datos: nivel vs. estructura

10 Glosario

Término Qué significa
AIC Medida para comparar modelos: premia el buen ajuste y castiga la complejidad. Más pequeño, mejor.
Binomial negativa Distribución para conteos que permite más variabilidad que la Poisson.
Coeficiente de variación Desviación estándar dividida por la media. Mide qué tan variables son los datos en relación con su promedio.
Devianza Qué tan lejos está el modelo de un ajuste perfecto. Más pequeña, mejor.
Deciles División de los datos ordenados en 10 grupos del mismo tamaño.
Diagnóstico Revisión de si el modelo cumple sus supuestos.
Error estándar (EE) Qué tan seguro está el modelo de un coeficiente. Pequeño = mucha seguridad.
Exposición Tiempo que una póliza estuvo vigente, en años.
Factor El exponencial de un coeficiente, e^β. Multiplica la predicción: 1,2 = +20 %; 0,8 = −20 %.
Frecuencia Número de siniestros por año de exposición.
Función de enlace El “traductor” entre la suma de variables y la escala de Y. En seguros suele ser el logaritmo.
Gaussiana inversa Distribución para montos positivos con una cola derecha aún más larga que la Gamma.
GLM Modelo Lineal Generalizado: regresión que permite elegir la distribución y la escala.
Intervalo de confianza Rango donde razonablemente está el verdadero valor de un parámetro.
Lift Gráfico que muestra si el modelo distingue a los clientes de bajo y alto costo.
Máxima verosimilitud Método para estimar coeficientes: elegir los valores que hacen más creíbles los datos observados.
Monitoreo Seguimiento del modelo después de implementado, para detectar si deja de funcionar.
Offset Variable que entra al modelo sin coeficiente, para ajustar por el tiempo de exposición.
p-valor Probabilidad de ver un resultado así si la variable no tuviera efecto. Menor a 0,05 suele indicar que sí aporta.
Prima comercial Precio que paga el cliente: prima pura más gastos y utilidad.
Prima pura Costo esperado de los siniestros de un cliente en un año: frecuencia × severidad.
Residual Diferencia entre lo observado y lo que predice el modelo.
Residual cuantílico Residual transformado para que, si el modelo es correcto, se comporte como una normal.
Severidad Costo promedio de cada siniestro.
Sobredispersión Cuando los datos varían más de lo que supone la distribución elegida.
Train / test Datos para entrenar el modelo y datos guardados para evaluarlo.
Tweedie Distribución para variables con muchos ceros y montos positivos, como el costo anual de una póliza.
VIF Indicador de si dos o más variables explicativas repiten la misma información.

11 Referencias

  • De Jong, P. y Heller, G. Z. (2008). Generalized Linear Models for Insurance Data. Cambridge University Press. (Fuente de la base dataCar.)
  • Dunn, P. K. y Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5(3), 236–244.
  • Dunn, P. K. y Smyth, G. K. (2018). Generalized Linear Models With Examples in R. Springer.
  • Goldburd, M., Khare, A., Tevet, D. y Guller, D. (2020). Generalized Linear Models for Insurance Rating (2.ª ed.). CAS Monograph Series n.º 5.
  • McCullagh, P. y Nelder, J. A. (1989). Generalized Linear Models (2.ª ed.). Chapman & Hall.
  • Nelder, J. A. y Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society, Series A, 135(3), 370–384.
  • Ohlsson, E. y Johansson, B. (2010). Non-Life Insurance Pricing with Generalized Linear Models. Springer.