Regresión Logística

Maestría en Ciencias Actuariales Aplicadas

1 Introducción

La regresión logística es una de las herramientas fundamentales para el análisis estadístico de fenómenos actuariales cuando la variable de interés es categórica.

En Ciencias Actuariales existen numerosos problemas en los que no se busca predecir directamente una cantidad monetaria, sino estimar la probabilidad de ocurrencia de un evento.

Algunos ejemplos son:

• determinar la probabilidad de que un asegurado presente al menos un siniestro;
• estimar la probabilidad de cancelación o no renovación de una póliza;
• identificar reclamaciones con mayor probabilidad de fraude;
• modelar la probabilidad de incumplimiento;
• clasificar asegurados de acuerdo con diferentes tipos de siniestro;
• estimar la probabilidad de ocurrencia de eventos de alta severidad;
• identificar segmentos de mayor riesgo para fines de suscripción o tarificación.

Supongamos que se desea estudiar si un asegurado presentará un siniestro durante el siguiente período de cobertura.

Definimos:

\[ Y_i = \begin{cases} 1, & \text{si el asegurado }i\text{ presenta al menos un siniestro},\\ 0, & \text{si el asegurado }i\text{ no presenta ningún siniestro}. \end{cases} \]

Nuestro objetivo consiste en estimar

\[ p_i=P(Y_i=1\mid X_i), \]

donde \(X_i\) representa las características conocidas del asegurado.

Por ejemplo:

\[ X_i = ( \text{edad}, \text{antigüedad}, \text{siniestros previos}, \text{zona}, \text{uso del vehículo}, \text{valor asegurado} ). \]

La regresión logística permite relacionar estas características con la probabilidad de ocurrencia del evento asegurado.

Idea actuarial central: la regresión logística no se limita a clasificar observaciones en 0 y 1. Su principal fortaleza consiste en producir una probabilidad individual de ocurrencia, que posteriormente puede utilizarse como insumo para decisiones actuariales.

1.1 Objetivos de aprendizaje

Al finalizar esta clase, el estudiante deberá ser capaz de:

1. Comprender la estructura matemática de la regresión logística.
2. Diferenciar probabilidad, odds y log-odds.
3. Comprender por qué una regresión lineal no resulta apropiada para una variable binaria.
4. Estimar un modelo logístico mediante glm().
5. Formular e interpretar contrastes de hipótesis.
6. Interpretar coeficientes logísticos.
7. Calcular e interpretar Odds Ratios.
8. Construir intervalos de confianza para los Odds Ratios.
9. Calcular probabilidades individuales de siniestro.
10. Evaluar el ajuste mediante pseudo \(R^2\).
11. Analizar la capacidad discriminante mediante ROC y AUC.
12. Evaluar la calibración mediante la prueba de Hosmer-Lemeshow.
13. Comprender la diferencia entre discriminación y calibración.
14. Ajustar una regresión logística multinomial.
15. Interpretar modelos multinomiales en un contexto de tipos de siniestro.

---

1.2 ¿Por qué no utilizar regresión lineal?

Supongamos que intentamos utilizar:

\[ Y_i= \beta_0+ \beta_1X_{1i}+ \cdots+ \beta_kX_{ki}+ \varepsilon_i, \]

donde

\[ Y_i\in\{0,1\}. \]

Este enfoque presenta varios inconvenientes.

Primero, una combinación lineal:

\[ \beta_0+ \beta_1X_1+ \cdots+ \beta_kX_k \]

puede tomar cualquier valor real.

Por lo tanto, una regresión lineal podría producir:

\[ \hat p=-0.10 \]

o

\[ \hat p=1.20, \]

aunque una probabilidad debe encontrarse siempre dentro del intervalo:

\[ 0\leq p\leq1. \]

Además, cuando la variable respuesta es binaria:

\[ Var(Y_i\mid X_i) = p_i(1-p_i), \]

por lo que la varianza depende de la propia probabilidad y no es constante.

La regresión logística resuelve estos problemas utilizando una transformación no lineal.

---

1.3 La función logística

Definimos el predictor lineal:

\[ \eta_i= \beta_0+ \beta_1X_{1i}+ \beta_2X_{2i} +\cdots+ \beta_kX_{ki}. \]

La función logística transforma este valor en una probabilidad:

\[ p_i= \frac{e^{\eta_i}} {1+e^{\eta_i}}. \]

Equivalentemente:

\[ p_i= \frac{1} {1+e^{-\eta_i}}. \]

Esta función garantiza:

\[ 0<p_i<1. \]

Observemos gráficamente la función logística.

eta_grafico <- seq(
  from = -8,
  to = 8,
  length.out = 500
)

datos_logistica <- tibble(
  eta = eta_grafico,
  probabilidad = plogis(eta_grafico)
)

ggplot(
  datos_logistica,
  aes(
    x = eta,
    y = probabilidad
  )
) +
  geom_line(
    linewidth = 1
  ) +
  geom_hline(
    yintercept = 0.5,
    linetype = "dashed"
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Función logística",
    subtitle = "Transformación del predictor lineal en una probabilidad",
    x = expression(eta),
    y = "Probabilidad"
  ) +
  theme_minimal()

1.3.1 Interpretación

Cuando:

\[ \eta=0, \]

se obtiene:

\[ p= \frac{e^0}{1+e^0} = \frac{1}{2} = 0.50. \]

Cuando \(\eta\) aumenta, la probabilidad se aproxima a 1.

Cuando \(\eta\) disminuye, la probabilidad se aproxima a 0.

Esta relación no es lineal.

Por ello, un mismo cambio en una variable explicativa no produce necesariamente el mismo cambio absoluto en probabilidad para todos los asegurados.

---

1.4 Probabilidad, odds y log-odds

La correcta interpretación de una regresión logística requiere comprender tres conceptos.

1.4.1 Probabilidad

Supongamos:

\[ p=0.20. \]

Esto significa que el evento tiene una probabilidad del 20%.

1.4.2 Odds

Los odds o momios se definen como:

\[ Odds= \frac{p}{1-p}. \]

Si:

\[ p=0.20, \]

entonces:

\[ Odds= \frac{0.20}{0.80} = 0.25. \]

Podemos interpretarlo como aproximadamente:

\[ 1:4, \]

es decir, un evento por cada cuatro no eventos.

Veamos varios ejemplos.

ejemplo_odds <- tibble(
  probabilidad = c(
    0.05,
    0.10,
    0.20,
    0.50,
    0.75,
    0.90
  )
) %>%
  mutate(
    odds = probabilidad /
      (1 - probabilidad)
  )

knitr::kable(
  ejemplo_odds,
  digits = 3,
  caption = "Relación entre probabilidad y odds"
)

Relación entre probabilidad y odds

probabilidad	odds
0.05	0.053
0.10	0.111
0.20	0.250
0.50	1.000
0.75	3.000
0.90	9.000

1.4.3 Log-odds o logit

La transformación logit se define como:

\[ \operatorname{logit}(p) = \log \left( \frac{p}{1-p} \right). \]

Por tanto, el modelo logístico puede escribirse como:

\[ \log \left( \frac{p_i}{1-p_i} \right) = \beta_0+ \beta_1X_{1i} +\cdots+ \beta_kX_{ki}. \]

Esta es la ecuación fundamental de la regresión logística binomial.

---

2 El modelo de regresión logística binomial

Desarrollaremos un caso completamente reproducible relacionado con una cartera ficticia de seguros de automóviles.

La variable de interés será:

\[ Y= \begin{cases} 1, & \text{si la póliza presenta al menos un siniestro},\\ 0, & \text{si no presenta siniestro}. \end{cases} \]

Utilizaremos como predictores:

• edad del asegurado;
• antigüedad de la póliza;
• número de siniestros previos;
• zona de riesgo;
• uso del vehículo;
• valor asegurado.

Todos los datos serán simulados dentro del documento.

No se requiere descargar ningún archivo externo.

---

2.1 Construcción reproducible de la cartera actuarial

Antes de estimar el modelo necesitamos una base de datos.

En una aplicación real esta información podría provenir de los sistemas de emisión, siniestros, cobranza o suscripción de una aseguradora.

En esta clase simularemos la información utilizando set.seed() para garantizar que cualquier estudiante obtenga los mismos resultados.

Supongamos que disponemos de 6,000 pólizas.

set.seed(10102026)

n <- 6000

seguros <- tibble(
  
  id_poliza = 1:n,
  
  edad = round(
    pmin(
      pmax(
        rnorm(
          n,
          mean = 42,
          sd = 13
        ),
        18
      ),
      80
    )
  ),
  
  antiguedad_poliza = round(
    pmin(
      rgamma(
        n,
        shape = 2.3,
        scale = 2.3
      ),
      20
    ),
    1
  ),
  
  siniestros_previos = pmin(
    rpois(
      n,
      lambda = 0.50
    ),
    5
  ),
  
  zona = sample(
    c(
      "Baja",
      "Media",
      "Alta"
    ),
    size = n,
    replace = TRUE,
    prob = c(
      0.38,
      0.42,
      0.20
    )
  ),
  
  uso_vehiculo = sample(
    c(
      "Particular",
      "Comercial"
    ),
    size = n,
    replace = TRUE,
    prob = c(
      0.82,
      0.18
    )
  ),
  
  valor_asegurado = round(
    exp(
      rnorm(
        n,
        mean = log(125000),
        sd = 0.45
      )
    )
  )
)

seguros <- seguros %>%
  mutate(
    
    zona = factor(
      zona,
      levels = c(
        "Baja",
        "Media",
        "Alta"
      )
    ),
    
    uso_vehiculo = factor(
      uso_vehiculo,
      levels = c(
        "Particular",
        "Comercial"
      )
    )
    
  )

head(seguros, 10)

2.1.1 Interpretación

Cada fila representa una póliza.

Por el momento únicamente disponemos de características del riesgo.

Necesitamos simular también si cada póliza presentó o no un siniestro.

Para generar este resultado construiremos una probabilidad subyacente de siniestro.

Supongamos que, en el proceso generador de los datos:

• la edad reduce moderadamente el riesgo;
• la antigüedad de la póliza tiene un pequeño efecto positivo;
• cada siniestro previo aumenta el riesgo;
• una zona de riesgo media presenta mayor riesgo que una zona baja;
• una zona alta presenta un incremento aún mayor;
• el uso comercial aumenta la probabilidad;
• un mayor valor asegurado tiene un efecto moderado.

El predictor lineal utilizado para simular los datos será:

\[ \eta_i= -3.15 -0.020(Edad_i-40) +0.035Antigüedad_i +0.65SiniestrosPrevios_i \]

\[ +0.40ZonaMedia_i +0.95ZonaAlta_i +0.60UsoComercial_i +0.0000025(ValorAsegurado_i-125000). \]

eta_real <- 
  -3.15 +
  (-0.020 * (seguros$edad - 40)) +
  (0.035 * seguros$antiguedad_poliza) +
  (0.65 * seguros$siniestros_previos) +
  (0.40 * (seguros$zona == "Media")) +
  (0.95 * (seguros$zona == "Alta")) +
  (0.60 * (seguros$uso_vehiculo == "Comercial")) +
  (0.0000025 * (seguros$valor_asegurado - 125000))

seguros <- seguros %>%
  mutate(
    
    probabilidad_real = plogis(
      eta_real
    ),
    
    siniestro = rbinom(
      n = n(),
      size = 1,
      prob = probabilidad_real
    )
    
  )

head(seguros, 10)

2.1.2 Interpretación

La variable probabilidad_real representa la probabilidad utilizada para simular los datos.

En una aplicación real esta probabilidad no sería observable.

Precisamente la finalidad de la regresión logística consiste en estimar:

\[ P(Y=1\mid X). \]

La variable siniestro contiene el resultado observado:

• 0: no ocurrió ningún siniestro;
• 1: ocurrió al menos un siniestro.

---

2.2 Exploración inicial del portafolio

Antes de estimar un modelo es indispensable realizar una exploración descriptiva.

summary(
  seguros %>%
    select(
      edad,
      antiguedad_poliza,
      siniestros_previos,
      zona,
      uso_vehiculo,
      valor_asegurado,
      siniestro
    )
)

##       edad      antiguedad_poliza siniestros_previos    zona     
##  Min.   :18.0   Min.   : 0.10     Min.   :0.000      Baja :2325  
##  1st Qu.:33.0   1st Qu.: 2.70     1st Qu.:0.000      Media:2490  
##  Median :42.0   Median : 4.60     Median :0.000      Alta :1185  
##  Mean   :42.2   Mean   : 5.29     Mean   :0.516                  
##  3rd Qu.:51.0   3rd Qu.: 7.10     3rd Qu.:1.000                  
##  Max.   :80.0   Max.   :20.00     Max.   :4.000                  
##      uso_vehiculo  valor_asegurado    siniestro    
##  Particular:4900   Min.   : 21656   Min.   :0.000  
##  Comercial :1100   1st Qu.: 92583   1st Qu.:0.000  
##                    Median :125264   Median :0.000  
##                    Mean   :138536   Mean   :0.118  
##                    3rd Qu.:170136   3rd Qu.:0.000  
##                    Max.   :599507   Max.   :1.000

Calculamos el número de pólizas y la frecuencia global del evento.

resumen_global <- seguros %>%
  summarise(
    polizas = n(),
    polizas_con_siniestro = sum(siniestro),
    polizas_sin_siniestro = sum(siniestro == 0),
    frecuencia_siniestro = mean(siniestro)
  )

knitr::kable(
  resumen_global,
  digits = 4,
  caption = "Resumen global de la cartera"
)

Resumen global de la cartera

polizas	polizas_con_siniestro	polizas_sin_siniestro	frecuencia_siniestro
6000	709	5291	0.1182

La frecuencia observada de pólizas con al menos un siniestro es:

\[ \hat p= 0.1182. \]

Expresada porcentualmente:

\[ 11.82%. \]

2.2.1 Interpretación actuarial

Si asignáramos una única probabilidad a toda la cartera, podríamos utilizar la frecuencia global:

\[ \hat p= \frac{ \text{Pólizas con siniestro} }{ \text{Total de pólizas} }. \]

Sin embargo, este procedimiento supondría que todos los asegurados tienen el mismo nivel de riesgo.

Una función esencial de la modelización actuarial consiste precisamente en reconocer la heterogeneidad de riesgos.

---

2.3 Frecuencia por zona

Analicemos la frecuencia observada de siniestros según la zona de riesgo.

frecuencia_zona <- seguros %>%
  group_by(zona) %>%
  summarise(
    polizas = n(),
    siniestros = sum(siniestro),
    frecuencia = mean(siniestro),
    .groups = "drop"
  )

knitr::kable(
  frecuencia_zona,
  digits = 4,
  caption = "Frecuencia observada de siniestro por zona"
)

Frecuencia observada de siniestro por zona

zona	polizas	siniestros	frecuencia
Baja	2325	205	0.0882
Media	2490	297	0.1193
Alta	1185	207	0.1747

Representamos los resultados.

ggplot(
  frecuencia_zona,
  aes(
    x = zona,
    y = frecuencia
  )
) +
  geom_col() +
  labs(
    title = "Frecuencia observada de siniestros por zona",
    x = "Zona de riesgo",
    y = "Frecuencia observada"
  ) +
  theme_minimal()

2.3.1 Interpretación

La comparación descriptiva permite observar si existen diferencias entre las zonas.

Sin embargo, todavía no podemos concluir que esas diferencias se deban exclusivamente a la zona.

Por ejemplo, los asegurados de una zona podrían tener también:

• mayor número de siniestros previos;
• vehículos de mayor valor;
• mayor proporción de uso comercial;
• edades diferentes.

La regresión logística permite estimar el efecto de una variable controlando simultáneamente las demás variables incluidas en el modelo.

---

2.4 Separación entre entrenamiento y validación

En aplicaciones predictivas es recomendable distinguir entre:

• datos utilizados para estimar el modelo;
• datos utilizados para evaluar su comportamiento.

Utilizaremos aproximadamente 70% para entrenamiento y 30% para validación.

set.seed(20261010)

indice_entrenamiento <- sample(
  seq_len(nrow(seguros)),
  size = floor(
    0.70 * nrow(seguros)
  )
)

entrenamiento <- seguros[
  indice_entrenamiento,
]

validacion <- seguros[
  -indice_entrenamiento,
]

tibble(
  conjunto = c(
    "Entrenamiento",
    "Validación"
  ),
  observaciones = c(
    nrow(entrenamiento),
    nrow(validacion)
  ),
  frecuencia_siniestro = c(
    mean(entrenamiento$siniestro),
    mean(validacion$siniestro)
  )
) %>%
  knitr::kable(
    digits = 4,
    caption = "Distribución de los datos entre entrenamiento y validación"
  )

Distribución de los datos entre entrenamiento y validación

conjunto	observaciones	frecuencia_siniestro
Entrenamiento	4200	0.1195
Validación	1800	0.1150

2.4.1 Interpretación

El modelo será estimado exclusivamente con el conjunto de entrenamiento.

Posteriormente utilizaremos el conjunto de validación para estudiar:

• discriminación;
• clasificación;
• calibración.

Esto permite evaluar el comportamiento del modelo en observaciones que no fueron utilizadas directamente para calcular sus coeficientes.

---

3 Estimación del modelo

En una regresión logística binomial suponemos:

\[ Y_i\sim Bernoulli(p_i). \]

Por lo tanto:

\[ P(Y_i=1)=p_i \]

y

\[ P(Y_i=0)=1-p_i. \]

La distribución Bernoulli puede escribirse como:

\[ P(Y_i=y_i) = p_i^{y_i} (1-p_i)^{1-y_i}. \]

La relación entre las variables explicativas y la probabilidad se plantea mediante:

\[ \log \left( \frac{p_i}{1-p_i} \right) = \beta_0+ \beta_1X_{1i} +\cdots+ \beta_kX_{ki}. \]

---

3.1 Estimación por máxima verosimilitud

En regresión logística no se utilizan mínimos cuadrados ordinarios.

Los coeficientes se estiman mediante máxima verosimilitud.

Para \(n\) observaciones independientes:

\[ L(\boldsymbol{\beta}) = \prod_{i=1}^{n} p_i^{y_i} (1-p_i)^{1-y_i}. \]

Es habitual trabajar con la log-verosimilitud:

\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[ y_i\log(p_i) + (1-y_i)\log(1-p_i) \right]. \]

Los valores:

\[ \hat\beta_0, \hat\beta_1, \ldots, \hat\beta_k \]

son aquellos que maximizan esta función.

---

3.2 Ajuste mediante glm()

En R, la regresión logística binomial puede estimarse mediante:

glm(
  formula,
  data,
  family = binomial(link = "logit")
)

Ajustaremos el modelo:

\[ Siniestro \sim Edad+ Antigüedad+ SiniestrosPrevios+ Zona+ Uso+ ValorAsegurado. \]

modelo_logit <- glm(
  
  siniestro ~
    edad +
    antiguedad_poliza +
    siniestros_previos +
    zona +
    uso_vehiculo +
    valor_asegurado,
  
  data = entrenamiento,
  
  family = binomial(
    link = "logit"
  )
  
)

summary(modelo_logit)

## 
## Call:
## glm(formula = siniestro ~ edad + antiguedad_poliza + siniestros_previos + 
##     zona + uso_vehiculo + valor_asegurado, family = binomial(link = "logit"), 
##     data = entrenamiento)
## 
## Coefficients:
##                           Estimate   Std. Error z value             Pr(>|z|)
## (Intercept)           -2.852235865  0.228142357  -12.50 < 0.0000000000000002
## edad                  -0.018359587  0.003932105   -4.67         0.0000030245
## antiguedad_poliza      0.034723063  0.013925420    2.49              0.01265
## siniestros_previos     0.686212590  0.059226296   11.59 < 0.0000000000000002
## zonaMedia              0.413562837  0.117250727    3.53              0.00042
## zonaAlta               0.804132547  0.130875381    6.14         0.0000000008
## uso_vehiculoComercial  0.532256421  0.117210794    4.54         0.0000055983
## valor_asegurado        0.000003429  0.000000692    4.95         0.0000007268
##                          
## (Intercept)           ***
## edad                  ***
## antiguedad_poliza     *  
## siniestros_previos    ***
## zonaMedia             ***
## zonaAlta              ***
## uso_vehiculoComercial ***
## valor_asegurado       ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3074.2  on 4199  degrees of freedom
## Residual deviance: 2846.0  on 4192  degrees of freedom
## AIC: 2862
## 
## Number of Fisher Scoring iterations: 5

3.2.1 ¿Qué información presenta summary()?

El resumen contiene, entre otros elementos:

• coeficientes estimados;
• errores estándar;
• estadísticos \(z\);
• p-valores;
• deviance nula;
• deviance residual;
• AIC.

La ecuación estimada tiene la forma:

\[ \widehat{\operatorname{logit}(p_i)} = \hat\beta_0+ \hat\beta_1Edad_i+ \hat\beta_2Antigüedad_i+ \hat\beta_3SiniestrosPrevios_i \]

\[ + \hat\beta_4ZonaMedia_i+ \hat\beta_5ZonaAlta_i+ \hat\beta_6UsoComercial_i+ \hat\beta_7ValorAsegurado_i. \]

Las categorías de referencia son:

• zona = Baja;
• uso_vehiculo = Particular.

---

3.3 Tabla profesional de coeficientes

Construimos una tabla más sencilla de leer.

tabla_coeficientes <- as.data.frame(
  summary(modelo_logit)$coefficients
) %>%
  rownames_to_column(
    "Variable"
  ) %>%
  rename(
    Coeficiente = Estimate,
    `Error estándar` = `Std. Error`,
    `Estadístico z` = `z value`,
    `p-valor` = `Pr(>|z|)`
  )

knitr::kable(
  tabla_coeficientes,
  digits = 4,
  caption = "Coeficientes estimados del modelo logístico"
)

Coeficientes estimados del modelo logístico

Variable	Coeficiente	Error estándar	Estadístico z	p-valor
(Intercept)	-2.8522	0.2281	-12.502	0.0000
edad	-0.0184	0.0039	-4.669	0.0000
antiguedad_poliza	0.0347	0.0139	2.494	0.0126
siniestros_previos	0.6862	0.0592	11.586	0.0000
zonaMedia	0.4136	0.1173	3.527	0.0004
zonaAlta	0.8041	0.1309	6.144	0.0000
uso_vehiculoComercial	0.5323	0.1172	4.541	0.0000
valor_asegurado	0.0000	0.0000	4.954	0.0000

3.3.1 Primera interpretación de los signos

Si:

\[ \hat\beta_j>0, \]

un incremento en \(X_j\) está asociado con un aumento en los log-odds de siniestro, manteniendo constantes las demás variables.

Si:

\[ \hat\beta_j<0, \]

un incremento en \(X_j\) está asociado con una reducción en los log-odds.

Esta interpretación todavía se encuentra en la escala logarítmica.

Más adelante transformaremos los coeficientes mediante:

\[ e^{\hat\beta_j} \]

para obtener Odds Ratios.

---

3.4 Obtención de probabilidades estimadas

Para cada póliza podemos calcular:

\[ \hat\eta_i = X_i^\top\hat\beta. \]

Luego:

\[ \hat p_i = \frac{ e^{\hat\eta_i} }{ 1+e^{\hat\eta_i} }. \]

Calculamos las probabilidades para el conjunto de entrenamiento.

entrenamiento <- entrenamiento %>%
  mutate(
    
    probabilidad_estimada = predict(
      modelo_logit,
      newdata = entrenamiento,
      type = "response"
    )
    
  )

entrenamiento %>%
  select(
    id_poliza,
    edad,
    siniestros_previos,
    zona,
    uso_vehiculo,
    siniestro,
    probabilidad_estimada
  ) %>%
  head(10) %>%
  knitr::kable(
    digits = 4,
    caption = "Ejemplo de probabilidades estimadas"
  )

Ejemplo de probabilidades estimadas

id_poliza	edad	siniestros_previos	zona	uso_vehiculo	siniestro	probabilidad_estimada
5208	42	0	Media	Particular	0	0.0599
65	35	0	Media	Particular	0	0.0826
979	28	0	Baja	Comercial	0	0.0743
1269	56	0	Alta	Particular	0	0.0712
859	57	1	Media	Particular	0	0.1099
872	45	1	Media	Comercial	0	0.1767
3961	48	0	Alta	Comercial	0	0.1567
3959	41	1	Media	Comercial	1	0.3459
3407	31	0	Baja	Particular	0	0.0492
4618	46	0	Media	Particular	0	0.0676

También calculamos probabilidades para el conjunto de validación.

validacion <- validacion %>%
  mutate(
    
    probabilidad_estimada = predict(
      modelo_logit,
      newdata = validacion,
      type = "response"
    )
    
  )

head(
  validacion %>%
    select(
      id_poliza,
      siniestro,
      probabilidad_estimada
    ),
  10
) %>%
  knitr::kable(
    digits = 4,
    caption = "Probabilidades estimadas en el conjunto de validación"
  )

Probabilidades estimadas en el conjunto de validación

id_poliza	siniestro	probabilidad_estimada
2	0	0.0773
4	0	0.0882
5	0	0.2073
6	0	0.1126
16	0	0.0374
27	0	0.1225
37	0	0.0664
42	0	0.1041
45	0	0.2074
48	0	0.0680

3.4.1 Interpretación actuarial

Si para una póliza obtenemos:

\[ \hat p_i=0.27, \]

el modelo estima una probabilidad de 27% de que dicha póliza presente al menos un siniestro durante el período considerado, condicionada a las variables incluidas.

Estas probabilidades pueden ser utilizadas como insumo para:

• segmentación;
• suscripción;
• tarificación;
• monitoreo de cartera;
• análisis de retención;
• prevención;
• priorización de inspecciones.

Debe tenerse presente que:

\[ P(\text{siniestro}) \]

no es equivalente a:

\[ E(\text{monto del siniestro}). \]

La frecuencia y la severidad representan dimensiones distintas del riesgo actuarial.

---

4 Contraste de hipótesis para el modelo estimado

Una vez estimado el modelo debemos determinar si las variables aportan evidencia estadística para explicar la probabilidad de siniestro.

Analizaremos:

1. contraste global mediante deviance;
2. contrastes individuales mediante estadísticos de Wald.

---

4.1 Deviance

La deviance es una medida relacionada con la log-verosimilitud.

En términos generales, una menor deviance representa un mejor ajuste relativo.

El summary() del modelo muestra:

• Null deviance;
• Residual deviance.

La deviance nula corresponde al modelo que únicamente contiene intercepto.

La deviance residual corresponde al modelo con los predictores incluidos.

---

4.2 Modelo nulo

Ajustamos:

\[ \operatorname{logit}(p_i) = \beta_0. \]

modelo_nulo <- glm(
  siniestro ~ 1,
  data = entrenamiento,
  family = binomial(
    link = "logit"
  )
)

summary(modelo_nulo)

## 
## Call:
## glm(formula = siniestro ~ 1, family = binomial(link = "logit"), 
##     data = entrenamiento)
## 
## Coefficients:
##             Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)  -1.9969     0.0476     -42 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3074.2  on 4199  degrees of freedom
## Residual deviance: 3074.2  on 4199  degrees of freedom
## AIC: 3076
## 
## Number of Fisher Scoring iterations: 4

El modelo nulo supone que todos los asegurados poseen la misma probabilidad estimada de siniestro.

Esto equivaldría a ignorar completamente la segmentación del riesgo.

---

4.3 Prueba global de razón de verosimilitudes

Queremos contrastar:

\[ H_0: \beta_1= \beta_2= \cdots= \beta_k=0 \]

contra:

\[ H_1: \text{al menos un }\beta_j\neq0. \]

En R:

comparacion_deviance <- anova(
  modelo_nulo,
  modelo_logit,
  test = "Chisq"
)

comparacion_deviance

Realizamos también el cálculo explícito.

deviance_nula <-
  modelo_logit$null.deviance

deviance_residual <-
  modelo_logit$deviance

reduccion_deviance <-
  deviance_nula -
  deviance_residual

gl_deviance <-
  modelo_logit$df.null -
  modelo_logit$df.residual

p_deviance <- pchisq(
  reduccion_deviance,
  df = gl_deviance,
  lower.tail = FALSE
)

tabla_deviance <- tibble(
  `Deviance nula` = deviance_nula,
  `Deviance residual` = deviance_residual,
  `Reducción de deviance` = reduccion_deviance,
  `Grados de libertad` = gl_deviance,
  `p-valor` = p_deviance
)

knitr::kable(
  tabla_deviance,
  digits = 4,
  caption = "Contraste global mediante reducción de deviance"
)

Contraste global mediante reducción de deviance

Deviance nula	Deviance residual	Reducción de deviance	Grados de libertad	p-valor
3074	2846	228.2	7	0

4.3.1 Interpretación

La reducción de deviance es:

\[ 228.21. \]

El p-valor es:

\[ p=< 0.0001. \]

Si:

\[ p<0.05, \]

rechazamos \(H_0\).

En ese caso concluimos que el conjunto de variables explicativas mejora significativamente el modelo respecto de uno que únicamente contiene intercepto.

4.3.2 Interpretación actuarial

En términos actuariales, un resultado significativo implica que existe evidencia de que las características del riesgo permiten diferenciar la probabilidad de siniestro entre pólizas.

Esto respalda estadísticamente la idea de segmentación de riesgos.

---

4.4 Contrastes de Wald

Para cada parámetro evaluamos:

\[ H_0:\beta_j=0 \]

contra:

\[ H_1:\beta_j\neq0. \]

El estadístico de Wald es:

\[ z_j= \frac{ \hat\beta_j }{ SE(\hat\beta_j) }. \]

Bajo \(H_0\):

\[ z_j\approx N(0,1). \]

R calcula automáticamente estos estadísticos.

tabla_wald <- tabla_coeficientes %>%
  mutate(
    `Significativo al 5%` =
      if_else(
        `p-valor` < 0.05,
        "Sí",
        "No"
      )
  )

knitr::kable(
  tabla_wald,
  digits = 4,
  caption = "Contrastes individuales de Wald"
)

Contrastes individuales de Wald

Variable	Coeficiente	Error estándar	Estadístico z	p-valor	Significativo al 5%
(Intercept)	-2.8522	0.2281	-12.502	0.0000	Sí
edad	-0.0184	0.0039	-4.669	0.0000	Sí
antiguedad_poliza	0.0347	0.0139	2.494	0.0126	Sí
siniestros_previos	0.6862	0.0592	11.586	0.0000	Sí
zonaMedia	0.4136	0.1173	3.527	0.0004	Sí
zonaAlta	0.8041	0.1309	6.144	0.0000	Sí
uso_vehiculoComercial	0.5323	0.1172	4.541	0.0000	Sí
valor_asegurado	0.0000	0.0000	4.954	0.0000	Sí

4.4.1 Interpretación

Si el p-valor de siniestros_previos es menor que 0.05, existe evidencia estadística de que el historial de siniestros está relacionado con la probabilidad futura de siniestro, manteniendo constantes las demás variables.

Sin embargo:

\[ \boxed{ \text{Significancia estadística} \neq \text{relevancia actuarial} } \]

Un efecto puede ser estadísticamente significativo y, al mismo tiempo, ser demasiado pequeño para modificar una decisión de tarificación o suscripción.

Por ello debemos estudiar conjuntamente:

• p-valores;
• magnitud del efecto;
• intervalos de confianza;
• Odds Ratios;
• coherencia técnica;
• estabilidad;
• desempeño predictivo.

---

5 Interpretación de los coeficientes de regresión

La interpretación correcta de los coeficientes es uno de los aspectos centrales del modelo logístico.

Recordemos:

\[ \log \left( \frac{p}{1-p} \right) = \beta_0+ \beta_1X_1+ \cdots+ \beta_kX_k. \]

Si \(X_j\) aumenta una unidad y mantenemos constantes las demás variables:

\[ \Delta\log(Odds) = \beta_j. \]

Al exponenciar:

\[ \frac{ Odds(X_j+1) }{ Odds(X_j) } = e^{\beta_j}. \]

Definimos:

\[ OR_j=e^{\beta_j}. \]

Esta cantidad recibe el nombre de Odds Ratio.

---

5.1 Cálculo de Odds Ratios

beta <- coef(
  modelo_logit
)

se_beta <- sqrt(
  diag(
    vcov(modelo_logit)
  )
)

tabla_or <- tibble(
  
  Variable = names(beta),
  
  Coeficiente = beta,
  
  `Odds Ratio` = exp(beta),
  
  `IC 95% inferior` =
    exp(
      beta -
        1.96 * se_beta
    ),
  
  `IC 95% superior` =
    exp(
      beta +
        1.96 * se_beta
    )
  
)

knitr::kable(
  tabla_or,
  digits = 4,
  caption = "Odds Ratios e intervalos de confianza del 95%"
)

Odds Ratios e intervalos de confianza del 95%

Variable	Coeficiente	Odds Ratio	IC 95% inferior	IC 95% superior
(Intercept)	-2.8522	0.0577	0.0369	0.0903
edad	-0.0184	0.9818	0.9743	0.9894
antiguedad_poliza	0.0347	1.0353	1.0075	1.0640
siniestros_previos	0.6862	1.9862	1.7685	2.2307
zonaMedia	0.4136	1.5122	1.2017	1.9029
zonaAlta	0.8041	2.2348	1.7291	2.8882
uso_vehiculoComercial	0.5323	1.7028	1.3533	2.1425
valor_asegurado	0.0000	1.0000	1.0000	1.0000

5.1.1 Regla general

Si:

\[ OR>1, \]

los odds del evento aumentan.

Si:

\[ OR<1, \]

los odds disminuyen.

Si:

\[ OR=1, \]

los odds no cambian.

---

5.2 Interpretación de siniestros previos

Extraemos el Odds Ratio.

or_siniestros_previos <- exp(
  coef(
    modelo_logit
  )[
    "siniestros_previos"
  ]
)

or_siniestros_previos

## siniestros_previos 
##              1.986

El valor estimado es:

\[ OR= 1.986. \]

El cambio porcentual en los odds es:

\[ 100(OR-1). \]

cambio_or_siniestros <-
  100 *
  (
    or_siniestros_previos -
      1
  )

cambio_or_siniestros

## siniestros_previos 
##              98.62

5.2.1 Interpretación actuarial

Manteniendo constantes:

• edad;
• antigüedad;
• zona;
• uso del vehículo;
• valor asegurado;

cada siniestro previo adicional multiplica los odds de presentar un nuevo siniestro por aproximadamente:

\[ 1.99. \]

Esto equivale a un cambio aproximado de:

\[ 98.6\% \]

en los odds.

Este resultado es consistente con el uso actuarial del historial de siniestros como indicador de experiencia del riesgo.

---

5.3 Interpretación de zona de riesgo

La categoría de referencia es:

\[ Zona=Baja. \]

Extraemos los Odds Ratios.

or_zona_media <- exp(
  coef(modelo_logit)[
    "zonaMedia"
  ]
)

or_zona_alta <- exp(
  coef(modelo_logit)[
    "zonaAlta"
  ]
)

tibble(
  Comparacion = c(
    "Zona Media vs. Zona Baja",
    "Zona Alta vs. Zona Baja"
  ),
  Odds_Ratio = c(
    or_zona_media,
    or_zona_alta
  )
) %>%
  knitr::kable(
    digits = 4,
    caption = "Odds Ratios para zona de riesgo"
  )

Odds Ratios para zona de riesgo

Comparacion	Odds_Ratio
Zona Media vs. Zona Baja	1.512
Zona Alta vs. Zona Baja	2.235

5.3.1 Interpretación

Para la zona alta:

\[ OR= 2.235. \]

Manteniendo constantes las demás variables, una póliza correspondiente a una zona de riesgo alta presenta odds de siniestro aproximadamente:

\[ 2.23 \]

veces los odds correspondientes a una póliza de zona baja.

---

5.4 Interpretación del uso comercial

La categoría de referencia es:

\[ Uso=Particular. \]

or_comercial <- exp(
  coef(modelo_logit)[
    "uso_vehiculoComercial"
  ]
)

or_comercial

## uso_vehiculoComercial 
##                 1.703

El Odds Ratio es:

\[ OR= 1.703. \]

5.4.1 Interpretación

Manteniendo constantes las demás variables, los vehículos de uso comercial presentan odds de siniestro aproximadamente:

\[ 1.70 \]

veces los correspondientes a vehículos de uso particular.

---

5.5 Odds Ratio no es Risk Ratio

Una confusión frecuente consiste en interpretar un Odds Ratio como si fuera directamente una razón de probabilidades.

Supongamos:

\[ p_1=0.20. \]

Entonces:

\[ Odds_1= \frac{0.20}{0.80} = 0.25. \]

Supongamos además:

\[ OR=2. \]

Los nuevos odds serían:

\[ Odds_2= 0.25(2) = 0.50. \]

La nueva probabilidad sería:

\[ p_2= \frac{0.50}{1+0.50} = 0.3333. \]

p_inicial <- 0.20

odds_inicial <-
  p_inicial /
  (1 - p_inicial)

OR_ejemplo <- 2

odds_nuevo <-
  odds_inicial *
  OR_ejemplo

p_nueva <-
  odds_nuevo /
  (1 + odds_nuevo)

tibble(
  `Probabilidad inicial` = p_inicial,
  `Odds iniciales` = odds_inicial,
  `Odds Ratio` = OR_ejemplo,
  `Nuevos odds` = odds_nuevo,
  `Nueva probabilidad` = p_nueva
) %>%
  knitr::kable(
    digits = 4,
    caption = "Diferencia entre Odds Ratio y cambio en probabilidad"
  )

Diferencia entre Odds Ratio y cambio en probabilidad

Probabilidad inicial	Odds iniciales	Odds Ratio	Nuevos odds	Nueva probabilidad
0.2	0.25	2	0.5	0.3333

5.5.1 Interpretación

Aunque los odds se duplicaron:

\[ OR=2, \]

la probabilidad no se duplicó.

Cambió de:

\[ 20\% \]

a aproximadamente:

\[ 33.3\%. \]

Por tanto:

\[ \boxed{ Odds\ Ratio \neq Risk\ Ratio } \]

---

5.6 Predicción para perfiles actuariales

Los Odds Ratios son útiles para estudiar efectos relativos.

Sin embargo, en la práctica también interesa responder preguntas como:

¿Cuál es la probabilidad estimada de siniestro de este asegurado?

Construiremos tres perfiles.

perfiles <- tibble(
  
  perfil = c(
    "Riesgo bajo",
    "Riesgo medio",
    "Riesgo alto"
  ),
  
  edad = c(
    55,
    40,
    25
  ),
  
  antiguedad_poliza = c(
    8,
    5,
    2
  ),
  
  siniestros_previos = c(
    0,
    1,
    3
  ),
  
  zona = factor(
    c(
      "Baja",
      "Media",
      "Alta"
    ),
    levels = levels(
      entrenamiento$zona
    )
  ),
  
  uso_vehiculo = factor(
    c(
      "Particular",
      "Particular",
      "Comercial"
    ),
    levels = levels(
      entrenamiento$uso_vehiculo
    )
  ),
  
  valor_asegurado = c(
    90000,
    140000,
    220000
  )
)

perfiles <- perfiles %>%
  mutate(
    
    probabilidad_siniestro = predict(
      modelo_logit,
      newdata = perfiles,
      type = "response"
    )
    
  )

knitr::kable(
  perfiles,
  digits = 4,
  caption = "Probabilidad estimada para diferentes perfiles actuariales"
)

Probabilidad estimada para diferentes perfiles actuariales

perfil	edad	antiguedad_poliza	siniestros_previos	zona	uso_vehiculo	valor_asegurado	probabilidad_siniestro
Riesgo bajo	55	8	0	Baja	Particular	90000	0.0364
Riesgo medio	40	5	1	Media	Particular	140000	0.1379
Riesgo alto	25	2	3	Alta	Comercial	220000	0.7125

Representamos las probabilidades.

ggplot(
  perfiles,
  aes(
    x = perfil,
    y = probabilidad_siniestro
  )
) +
  geom_col() +
  labs(
    title = "Probabilidad estimada por perfil actuarial",
    x = "Perfil",
    y = "Probabilidad estimada de siniestro"
  ) +
  theme_minimal()

5.6.1 Interpretación

La regresión logística permite transformar características observables en probabilidades individuales.

Este enfoque puede contribuir a:

• clasificación tarifaria;
• selección de riesgos;
• evaluación de políticas de suscripción;
• segmentación;
• monitoreo preventivo.

Debe enfatizarse que una probabilidad estimada es una medida de riesgo, no una decisión automática.

---

6 Evaluando el ajuste del modelo

En regresión lineal es habitual utilizar:

\[ R^2. \]

En regresión logística no existe una medida única con exactamente la misma interpretación.

Por esta razón se utilizan diferentes medidas denominadas pseudo \(R^2\).

---

6.1 Pseudo \(R^2\) de McFadden

El pseudo \(R^2\) de McFadden se define como:

\[ R^2_{McFadden} = 1- \frac{ \ell(\text{modelo completo}) }{ \ell(\text{modelo nulo}) }. \]

Calculamos:

loglik_completo <- as.numeric(
  logLik(
    modelo_logit
  )
)

loglik_nulo <- as.numeric(
  logLik(
    modelo_nulo
  )
)

r2_mcfadden <-
  1 -
  (
    loglik_completo /
      loglik_nulo
  )

tabla_mcfadden <- tibble(
  `LogLik modelo nulo` = loglik_nulo,
  `LogLik modelo completo` = loglik_completo,
  `Pseudo R2 McFadden` = r2_mcfadden
)

knitr::kable(
  tabla_mcfadden,
  digits = 4,
  caption = "Pseudo R-cuadrado de McFadden"
)

Pseudo R-cuadrado de McFadden

LogLik modelo nulo	LogLik modelo completo	Pseudo R2 McFadden
-1537	-1423	0.0742

El resultado es:

\[ R^2_{McFadden} = 0.0742. \]

6.1.1 Interpretación

Este valor indica cuánto mejora el modelo completo respecto del modelo que contiene únicamente intercepto.

No debe interpretarse como:

El modelo explica 7.42% de la variabilidad.

Esa interpretación no es válida.

El pseudo \(R^2\) de McFadden no es directamente comparable con el \(R^2\) de una regresión lineal.

---

6.2 Criterio de información AIC

El criterio de información de Akaike se define como:

\[ AIC= -2\ell+2k, \]

donde:

• \(\ell\) es la log-verosimilitud;
• \(k\) es el número de parámetros.

tabla_aic <- tibble(
  Modelo = c(
    "Modelo nulo",
    "Modelo logístico completo"
  ),
  AIC = c(
    AIC(modelo_nulo),
    AIC(modelo_logit)
  )
)

knitr::kable(
  tabla_aic,
  digits = 2,
  caption = "Comparación del AIC"
)

Comparación del AIC

Modelo	AIC
Modelo nulo	3076
Modelo logístico completo	2862

6.2.1 Interpretación

Cuando se comparan modelos:

• estimados con los mismos datos;
• para la misma variable respuesta;

un menor AIC indica un mejor compromiso entre:

• ajuste;
• complejidad.

Sin embargo, el AIC no mide directamente:

• discriminación;
• calibración;
• utilidad económica.

---

7 Capacidad discriminante del modelo

Un modelo puede ajustar adecuadamente los datos y, aun así, no separar correctamente a los asegurados con mayor y menor riesgo.

La discriminación estudia la capacidad del modelo para ordenar las observaciones de acuerdo con su propensión al evento.

Una pregunta actuarial sería:

Si seleccionamos aleatoriamente una póliza con siniestro y una póliza sin siniestro, ¿el modelo tiende a asignar una probabilidad mayor a la póliza que realmente tuvo el siniestro?

---

7.1 Clasificación mediante un umbral

Podemos transformar una probabilidad en una clasificación utilizando un punto de corte \(c\):

\[ \hat Y_i= \begin{cases} 1, & \hat p_i\geq c,\\ 0, & \hat p_i<c. \end{cases} \]

Por ejemplo:

\[ c=0.20. \]

Sin embargo, el punto de corte no altera las probabilidades producidas por el modelo; únicamente altera la regla utilizada para convertirlas en categorías.

---

7.2 Sensibilidad

La sensibilidad se define como:

\[ Sensibilidad = \frac{VP}{VP+FN}. \]

También:

\[ Sensibilidad = P(\hat Y=1\mid Y=1). \]

Responde:

De todos los asegurados que realmente presentaron un siniestro, ¿qué proporción fue identificada como positiva?

---

7.3 Especificidad

La especificidad es:

\[ Especificidad = \frac{VN}{VN+FP}. \]

También:

\[ Especificidad = P(\hat Y=0\mid Y=0). \]

Responde:

De todos los asegurados sin siniestro, ¿qué proporción fue identificada correctamente como negativa?

---

7.4 Curva ROC

La curva ROC representa:

\[ Sensibilidad \]

frente a:

\[ 1-Especificidad \]

para diferentes puntos de corte.

Evaluaremos la curva ROC utilizando el conjunto de validación.

roc_validacion <- roc(
  
  response =
    validacion$siniestro,
  
  predictor =
    validacion$probabilidad_estimada,
  
  levels = c(
    0,
    1
  ),
  
  direction = "<",
  
  quiet = TRUE
  
)

plot(
  roc_validacion,
  main = "Curva ROC - Modelo logístico",
  legacy.axes = TRUE,
  print.auc = TRUE,
  print.auc.cex = 1.1,
  lwd = 2
)

abline(
  a = 0,
  b = 1,
  lty = 2
)

---

7.5 Área bajo la curva: AUC

Calculamos:

auc_validacion <- as.numeric(
  auc(
    roc_validacion
  )
)

auc_validacion

## [1] 0.7149

El valor obtenido es:

\[ AUC= 0.715. \]

7.5.1 Interpretación probabilística

El AUC puede interpretarse aproximadamente como la probabilidad de que el modelo asigne un score de riesgo mayor a una observación con evento que a una observación sin evento seleccionadas aleatoriamente.

Por ejemplo, si:

\[ AUC=0.75, \]

podríamos interpretar que, en aproximadamente 75% de las parejas formadas por un siniestro y un no siniestro, el modelo ordena correctamente ambos riesgos.

7.5.2 Referencia orientativa

AUC	Interpretación aproximada
0.50	Sin discriminación
0.60 - 0.70	Discriminación limitada
0.70 - 0.80	Discriminación aceptable
0.80 - 0.90	Buena discriminación
> 0.90	Discriminación muy alta

Estas categorías son únicamente orientativas.

No existe un AUC universalmente aceptable para todos los problemas actuariales.

---

7.6 Punto de corte según índice de Youden

El índice de Youden se define como:

\[ J= Sensibilidad+ Especificidad- 1. \]

El punto que maximiza \(J\) puede utilizarse como referencia estadística.

punto_youden <- coords(
  roc_validacion,
  x = "best",
  best.method = "youden",
  ret = c(
    "threshold",
    "sensitivity",
    "specificity"
  )
)

punto_youden

7.6.1 Interpretación actuarial

El punto de corte de Youden no necesariamente representa la decisión económicamente óptima.

Por ejemplo, en detección de fraude:

• un falso positivo genera costos de investigación;
• un falso negativo puede generar una pérdida considerable.

Podríamos conceptualizar una función de costo como:

\[ Costo(c) = C_{FP}\times FP(c) + C_{FN}\times FN(c). \]

Por tanto, la selección del punto de corte debe considerar:

• costo del falso positivo;
• costo del falso negativo;
• capacidad operativa;
• apetito de riesgo;
• objetivos del modelo.

---

7.7 Matriz de confusión

Utilizaremos un umbral ilustrativo de 0.20.

umbral <- 0.20

validacion <- validacion %>%
  mutate(
    
    clasificacion = if_else(
      probabilidad_estimada >= umbral,
      1L,
      0L
    )
    
  )

matriz_confusion <- table(
  Observado = validacion$siniestro,
  Predicho = validacion$clasificacion
)

matriz_confusion

##          Predicho
## Observado    0    1
##         0 1440  153
##         1  144   63

Extraemos:

• verdaderos positivos;
• verdaderos negativos;
• falsos positivos;
• falsos negativos.

VP <- matriz_confusion[
  "1",
  "1"
]

VN <- matriz_confusion[
  "0",
  "0"
]

FP <- matriz_confusion[
  "0",
  "1"
]

FN <- matriz_confusion[
  "1",
  "0"
]

sensibilidad <-
  VP /
  (VP + FN)

especificidad <-
  VN /
  (VN + FP)

precision <-
  VP /
  (VP + FP)

exactitud <-
  (VP + VN) /
  sum(matriz_confusion)

tabla_metricas <- tibble(
  Metrica = c(
    "Sensibilidad",
    "Especificidad",
    "Precisión",
    "Exactitud"
  ),
  Valor = c(
    sensibilidad,
    especificidad,
    precision,
    exactitud
  )
)

knitr::kable(
  tabla_metricas,
  digits = 4,
  caption = paste(
    "Métricas de clasificación con umbral =",
    umbral
  )
)

Métricas de clasificación con umbral = 0.2

Metrica	Valor
Sensibilidad	0.3043
Especificidad	0.9040
Precisión	0.2917
Exactitud	0.8350

7.7.1 Interpretación

Para este punto de corte obtenemos aproximadamente:

• sensibilidad: 30.43%;
• especificidad: 90.40%;
• precisión: 29.17%;
• exactitud: 83.50%.

La exactitud debe analizarse con precaución cuando existe desbalance de clases.

Por ejemplo, si solamente 10% de las pólizas presentan siniestro, un modelo que predijera siempre:

No siniestro

tendría aproximadamente 90% de exactitud, aunque no identificaría ningún siniestro.

---

8 Calibración del modelo

La discriminación y la calibración responden preguntas diferentes.

8.1 Discriminación

Pregunta:

¿El modelo ordena correctamente los riesgos?

Una herramienta principal es:

\[ AUC. \]

8.2 Calibración

Pregunta:

¿Las probabilidades estimadas coinciden razonablemente con las frecuencias realmente observadas?

Por ejemplo, si un conjunto de asegurados recibe:

\[ \hat p\approx0.20, \]

esperaríamos observar aproximadamente 20% de eventos en ese grupo si las probabilidades están bien calibradas.

---

8.3 Calibración por deciles

Agruparemos las pólizas del conjunto de validación en diez grupos ordenados por riesgo estimado.

calibracion_deciles <- validacion %>%
  mutate(
    
    decil_riesgo = ntile(
      probabilidad_estimada,
      10
    )
    
  ) %>%
  group_by(
    decil_riesgo
  ) %>%
  summarise(
    
    polizas = n(),
    
    probabilidad_media =
      mean(
        probabilidad_estimada
      ),
    
    frecuencia_observada =
      mean(
        siniestro
      ),
    
    siniestros_observados =
      sum(
        siniestro
      ),
    
    siniestros_esperados =
      sum(
        probabilidad_estimada
      ),
    
    .groups = "drop"
    
  )

knitr::kable(
  calibracion_deciles,
  digits = 4,
  caption = "Calibración por deciles de riesgo"
)

Calibración por deciles de riesgo

decil_riesgo	polizas	probabilidad_media	frecuencia_observada	siniestros_observados	siniestros_esperados
1	180	0.0373	0.0278	5	6.717
2	180	0.0509	0.0389	7	9.166
3	180	0.0620	0.0389	7	11.158
4	180	0.0744	0.0611	11	13.397
5	180	0.0878	0.0722	13	15.813
6	180	0.1025	0.1444	26	18.451
7	180	0.1210	0.1444	26	21.784
8	180	0.1445	0.1278	23	26.016
9	180	0.1854	0.1778	32	33.364
10	180	0.3062	0.3167	57	55.118

---

8.4 Gráfico de calibración

ggplot(
  calibracion_deciles,
  aes(
    x = probabilidad_media,
    y = frecuencia_observada
  )
) +
  geom_point(
    size = 3
  ) +
  geom_line(
    linewidth = 0.8
  ) +
  geom_abline(
    intercept = 0,
    slope = 1,
    linetype = "dashed"
  ) +
  labs(
    title = "Gráfico de calibración",
    subtitle = "Probabilidad estimada frente a frecuencia observada por deciles",
    x = "Probabilidad media estimada",
    y = "Frecuencia observada"
  ) +
  theme_minimal()

8.4.1 Interpretación

La línea diagonal representa:

\[ Probabilidad\ estimada = Frecuencia\ observada. \]

Un modelo perfectamente calibrado tendría puntos muy próximos a esta diagonal.

Si los puntos se encuentran sistemáticamente por encima de la diagonal, el modelo tiende a subestimar el riesgo.

Si se encuentran sistemáticamente por debajo, tiende a sobreestimar el riesgo.

---

8.5 Test de Hosmer-Lemeshow

La prueba de Hosmer-Lemeshow compara:

• número observado de eventos;
• número esperado de eventos;

dentro de grupos definidos según la probabilidad estimada.

Podemos representar conceptualmente las hipótesis como:

\[ H_0: \text{no existe evidencia de discrepancias importantes entre observados y esperados} \]

contra:

\[ H_1: \text{existen discrepancias entre observados y esperados}. \]

Aplicamos la prueba sobre el conjunto de validación.

hl <- hoslem.test(
  
  x = validacion$siniestro,
  
  y = validacion$probabilidad_estimada,
  
  g = 10
  
)

hl

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  validacion$siniestro, validacion$probabilidad_estimada
## X-squared = 8.6, df = 8, p-value = 0.4

Extraemos el p-valor.

p_hl <- hl$p.value

p_hl

## [1] 0.3773

El p-valor obtenido es:

\[ p= 0.3773. \]

8.5.1 Interpretación

Si:

\[ p>0.05, \]

no rechazamos \(H_0\).

Esto significa que la prueba no proporciona evidencia estadísticamente significativa de falta de ajuste.

Si:

\[ p<0.05, \]

rechazamos \(H_0\).

Esto sugiere discrepancias estadísticamente detectables entre los resultados observados y esperados.

---

8.6 Precauciones con Hosmer-Lemeshow

La prueba de Hosmer-Lemeshow no debe utilizarse como una regla mecánica.

Su resultado puede depender de:

• tamaño de la muestra;
• número de grupos;
• distribución de las probabilidades;
• especificación del modelo.

En muestras grandes, desviaciones pequeñas pueden resultar estadísticamente significativas.

Por tanto, es recomendable complementar la prueba con:

1. gráficos de calibración;
2. comparación de observado versus esperado;
3. análisis por segmentos;
4. validación fuera de muestra;
5. validación temporal.

---

8.7 Razón observado/esperado

En seguros resulta intuitivo comparar:

\[ O/E= \frac{ Eventos\ observados }{ Eventos\ esperados }. \]

Calculamos:

observados <-
  sum(
    validacion$siniestro
  )

esperados <-
  sum(
    validacion$probabilidad_estimada
  )

razon_oe <-
  observados /
  esperados

tabla_oe <- tibble(
  `Siniestros observados` = observados,
  `Siniestros esperados` = esperados,
  `Razón O/E` = razon_oe
)

knitr::kable(
  tabla_oe,
  digits = 4,
  caption = "Comparación entre siniestros observados y esperados"
)

Comparación entre siniestros observados y esperados

Siniestros observados	Siniestros esperados	Razón O/E
207	211	0.9811

8.7.1 Interpretación

Si:

\[ O/E\approx1, \]

el número total de eventos pronosticados es similar al observado.

Si:

\[ O/E>1, \]

el modelo está subestimando el número total de eventos.

Si:

\[ O/E<1, \]

el modelo está sobreestimando el número total de eventos.

En nuestro ejemplo:

\[ O/E= 0.981. \]

---

8.8 Discriminación no es calibración

Esta diferencia es esencial.

Un modelo puede tener:

\[ AUC=0.85 \]

y, al mismo tiempo, producir probabilidades mal calibradas.

Supongamos que distingue perfectamente quién tiene mayor riesgo, pero asigna:

\[ 0.60 \]

a grupos cuya frecuencia real es:

\[ 0.30. \]

El ranking puede ser adecuado, pero las probabilidades no.

En una aplicación actuarial, esta diferencia tiene implicaciones importantes.

Para tarificación o estimación de pérdida esperada, necesitamos no solamente:

\[ \text{ordenar los riesgos}, \]

sino también:

\[ \text{estimar adecuadamente la magnitud de la probabilidad}. \]

---

9 Regresión logística multinomial

Hasta ahora hemos estudiado una respuesta binaria:

\[ Y\in\{0,1\}. \]

Sin embargo, muchos problemas actuariales presentan más de dos categorías.

Supongamos que una póliza puede finalizar el período en alguno de los siguientes estados:

\[ Y= \begin{cases} 0, & \text{Sin siniestro},\\ 1, & \text{Daños materiales},\\ 2, & \text{Lesiones},\\ 3, & \text{Robo}. \end{cases} \]

Estas categorías representan resultados mutuamente excluyentes.

Para este problema podemos utilizar una regresión logística multinomial.

---

9.1 Fundamento del modelo multinomial

Supongamos:

\[ Y_i\in \{ 0,1,\ldots,J-1 \}. \]

Seleccionamos una categoría como referencia.

En nuestro ejemplo utilizaremos:

\[ Y=0: \text{Sin siniestro}. \]

Para cada categoría \(j\neq0\):

\[ \log \left( \frac{ P(Y_i=j) }{ P(Y_i=0) } \right) = \beta_{j0} + \beta_{j1}X_{1i} +\cdots+ \beta_{jk}X_{ki}. \]

Por ejemplo:

\[ \log \left( \frac{ P(\text{Robo}) }{ P(\text{Sin siniestro}) } \right) = \beta_{R0} + \beta_{R1}Edad + \beta_{R2}SiniestrosPrevios +\cdots \]

La regresión multinomial estima una ecuación diferente para cada categoría distinta de la categoría base.

---

9.2 Simulación de datos multinomiales

Construiremos otra cartera ficticia.

set.seed(15102026)

n_multi <- 7000

seguros_multi <- tibble(
  
  id_poliza = 1:n_multi,
  
  edad = round(
    pmin(
      pmax(
        rnorm(
          n_multi,
          mean = 42,
          sd = 13
        ),
        18
      ),
      80
    )
  ),
  
  siniestros_previos = pmin(
    rpois(
      n_multi,
      lambda = 0.50
    ),
    5
  ),
  
  zona = sample(
    c(
      "Baja",
      "Media",
      "Alta"
    ),
    size = n_multi,
    replace = TRUE,
    prob = c(
      0.40,
      0.40,
      0.20
    )
  ),
  
  uso_vehiculo = sample(
    c(
      "Particular",
      "Comercial"
    ),
    size = n_multi,
    replace = TRUE,
    prob = c(
      0.83,
      0.17
    )
  ),
  
  valor_asegurado = round(
    exp(
      rnorm(
        n_multi,
        mean = log(125000),
        sd = 0.45
      )
    )
  )
)

seguros_multi <- seguros_multi %>%
  mutate(
    
    zona = factor(
      zona,
      levels = c(
        "Baja",
        "Media",
        "Alta"
      )
    ),
    
    uso_vehiculo = factor(
      uso_vehiculo,
      levels = c(
        "Particular",
        "Comercial"
      )
    )
    
  )

Definimos un predictor lineal diferente para cada tipo de siniestro.

eta_danos <-
  -2.40 +
  (-0.012 * (seguros_multi$edad - 40)) +
  (0.55 * seguros_multi$siniestros_previos) +
  (0.30 * (seguros_multi$zona == "Media")) +
  (0.55 * (seguros_multi$zona == "Alta")) +
  (0.35 * (seguros_multi$uso_vehiculo == "Comercial"))

eta_lesiones <-
  -3.70 +
  (-0.025 * (seguros_multi$edad - 40)) +
  (0.45 * seguros_multi$siniestros_previos) +
  (0.35 * (seguros_multi$zona == "Media")) +
  (0.80 * (seguros_multi$zona == "Alta")) +
  (0.65 * (seguros_multi$uso_vehiculo == "Comercial"))

eta_robo <-
  -3.20 +
  (-0.005 * (seguros_multi$edad - 40)) +
  (0.25 * seguros_multi$siniestros_previos) +
  (0.55 * (seguros_multi$zona == "Media")) +
  (1.35 * (seguros_multi$zona == "Alta")) +
  (0.20 * (seguros_multi$uso_vehiculo == "Comercial")) +
  (0.000003 * (seguros_multi$valor_asegurado - 125000))

---

9.3 Probabilidades multinomiales

La probabilidad de la categoría de referencia es:

\[ P(Y=0) = \frac{ 1 }{ 1+ e^{\eta_1} + e^{\eta_2} + e^{\eta_3} }. \]

Para una categoría \(j\):

\[ P(Y=j) = \frac{ e^{\eta_j} }{ 1+ \sum_{k=1}^{J-1}e^{\eta_k} }. \]

Calculamos:

denominador <-
  1 +
  exp(eta_danos) +
  exp(eta_lesiones) +
  exp(eta_robo)

p_sin <-
  1 /
  denominador

p_danos <-
  exp(eta_danos) /
  denominador

p_lesiones <-
  exp(eta_lesiones) /
  denominador

p_robo <-
  exp(eta_robo) /
  denominador

Verificamos que las probabilidades sumen 1.

suma_probabilidades <-
  p_sin +
  p_danos +
  p_lesiones +
  p_robo

summary(
  suma_probabilidades
)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1

Las pequeñas diferencias respecto de 1, si aparecen, se deben únicamente a precisión numérica.

---

9.4 Generación del tipo de siniestro

set.seed(16102026)

categorias <- c(
  "Sin siniestro",
  "Danos materiales",
  "Lesiones",
  "Robo"
)

tipo_siniestro <- map_chr(
  
  seq_len(
    n_multi
  ),
  
  function(i) {
    
    sample(
      
      categorias,
      
      size = 1,
      
      prob = c(
        p_sin[i],
        p_danos[i],
        p_lesiones[i],
        p_robo[i]
      )
      
    )
    
  }
  
)

seguros_multi <- seguros_multi %>%
  mutate(
    
    tipo_siniestro = factor(
      tipo_siniestro,
      levels = c(
        "Sin siniestro",
        "Danos materiales",
        "Lesiones",
        "Robo"
      )
    )
    
  )

head(seguros_multi, 10)

---

9.5 Distribución de los tipos de siniestro

distribucion_multi <- seguros_multi %>%
  count(
    tipo_siniestro
  ) %>%
  mutate(
    proporcion =
      n /
      sum(n)
  )

knitr::kable(
  distribucion_multi,
  digits = 4,
  caption = "Distribución de los tipos de resultado"
)

Distribución de los tipos de resultado

tipo_siniestro	n	proporcion
Sin siniestro	5312	0.7589
Danos materiales	901	0.1287
Lesiones	273	0.0390
Robo	514	0.0734

9.5.1 Interpretación

La regresión multinomial resulta apropiada porque la variable respuesta ya no contiene únicamente:

\[ 0/1. \]

Ahora cada póliza pertenece a una de varias categorías.

Además, diferentes variables pueden afectar de forma distinta:

• daños materiales;
• lesiones;
• robo.

---

9.6 Estimación mediante multinom()

La función multinom() pertenece al paquete nnet.

modelo_multinomial <- multinom(
  
  tipo_siniestro ~
    edad +
    siniestros_previos +
    zona +
    uso_vehiculo +
    valor_asegurado,
  
  data = seguros_multi,
  
  trace = FALSE
  
)

summary(
  modelo_multinomial
)

## Call:
## multinom(formula = tipo_siniestro ~ edad + siniestros_previos + 
##     zona + uso_vehiculo + valor_asegurado, data = seguros_multi, 
##     trace = FALSE)
## 
## Coefficients:
##                  (Intercept)      edad siniestros_previos zonaMedia zonaAlta
## Danos materiales      -2.053 -0.012009             0.5191    0.3065   0.3921
## Lesiones              -2.726 -0.028483             0.4890    0.4546   0.5362
## Robo                  -3.248 -0.007973             0.2646    0.6333   1.3572
##                  uso_vehiculoComercial valor_asegurado
## Danos materiales                0.3418     0.000001540
## Lesiones                        0.8118     0.000001114
## Robo                           -0.0757     0.000003492
## 
## Std. Errors:
##                  (Intercept)       edad siniestros_previos    zonaMedia
## Danos materiales 0.000026339 0.00162232        0.000019356 0.0000121911
## Lesiones         0.000001553 0.00009142        0.000001311 0.0000007821
## Robo             0.000006760 0.00041134        0.000004420 0.0000029417
##                      zonaAlta uso_vehiculoComercial valor_asegurado
## Danos materiales 0.0000057477          0.0000055475    0.0000004533
## Lesiones         0.0000003269          0.0000004907    0.0000004074
## Robo             0.0000024598          0.0000010506    0.0000002937
## 
## Residual Deviance: 10669 
## AIC: 10711

9.6.1 Interpretación

Debido a que la primera categoría es:

\[ \text{Sin siniestro}, \]

esta actúa como categoría de referencia.

El modelo estima:

\[ \log \left( \frac{ P(\text{Daños materiales}) }{ P(\text{Sin siniestro}) } \right), \]

\[ \log \left( \frac{ P(\text{Lesiones}) }{ P(\text{Sin siniestro}) } \right), \]

y:

\[ \log \left( \frac{ P(\text{Robo}) }{ P(\text{Sin siniestro}) } \right). \]

Por tanto, un mismo predictor puede presentar coeficientes distintos según el tipo de siniestro.

---

9.7 Pruebas de Wald en regresión multinomial

multinom() devuelve:

• coeficientes;
• errores estándar.

Podemos construir:

\[ z= \frac{ \hat\beta }{ SE(\hat\beta) } \]

y calcular p-valores aproximados.

resumen_multi <- summary(
  modelo_multinomial
)

coef_multi <-
  resumen_multi$coefficients

se_multi <-
  resumen_multi$standard.errors

z_multi <-
  coef_multi /
  se_multi

p_multi <-
  2 *
  (
    1 -
      pnorm(
        abs(
          z_multi
        )
      )
  )

Organizamos los resultados.

tabla_multi <- map_dfr(
  
  rownames(
    coef_multi
  ),
  
  function(categoria) {
    
    tibble(
      
      Categoria = categoria,
      
      Variable =
        colnames(
          coef_multi
        ),
      
      Coeficiente =
        coef_multi[
          categoria,
        ],
      
      `Error estándar` =
        se_multi[
          categoria,
        ],
      
      `Estadístico z` =
        z_multi[
          categoria,
        ],
      
      `p-valor` =
        p_multi[
          categoria,
        ]
      
    )
    
  }
  
)

knitr::kable(
  tabla_multi,
  digits = 4,
  caption = "Coeficientes y pruebas de Wald del modelo multinomial"
)

Coeficientes y pruebas de Wald del modelo multinomial

Categoria	Variable	Coeficiente	Error estándar	Estadístico z	p-valor
Danos materiales	(Intercept)	-2.0531	0.0000	-77947.486	0.0000
Danos materiales	edad	-0.0120	0.0016	-7.403	0.0000
Danos materiales	siniestros_previos	0.5191	0.0000	26817.937	0.0000
Danos materiales	zonaMedia	0.3065	0.0000	25139.210	0.0000
Danos materiales	zonaAlta	0.3921	0.0000	68219.450	0.0000
Danos materiales	uso_vehiculoComercial	0.3418	0.0000	61610.926	0.0000
Danos materiales	valor_asegurado	0.0000	0.0000	3.396	0.0007
Lesiones	(Intercept)	-2.7258	0.0000	-1754664.217	0.0000
Lesiones	edad	-0.0285	0.0001	-311.568	0.0000
Lesiones	siniestros_previos	0.4890	0.0000	373055.226	0.0000
Lesiones	zonaMedia	0.4546	0.0000	581276.075	0.0000
Lesiones	zonaAlta	0.5362	0.0000	1640203.725	0.0000
Lesiones	uso_vehiculoComercial	0.8118	0.0000	1654613.165	0.0000
Lesiones	valor_asegurado	0.0000	0.0000	2.734	0.0063
Robo	(Intercept)	-3.2477	0.0000	-480404.196	0.0000
Robo	edad	-0.0080	0.0004	-19.384	0.0000
Robo	siniestros_previos	0.2646	0.0000	59848.304	0.0000
Robo	zonaMedia	0.6333	0.0000	215283.553	0.0000
Robo	zonaAlta	1.3572	0.0000	551756.998	0.0000
Robo	uso_vehiculoComercial	-0.0757	0.0000	-72060.303	0.0000
Robo	valor_asegurado	0.0000	0.0000	11.890	0.0000

9.7.1 Interpretación

Cada coeficiente debe interpretarse con respecto a:

\[ \text{Sin siniestro}. \]

Por ejemplo, si el coeficiente de zonaAlta en la ecuación de Robo es positivo, una zona alta incrementa:

\[ \log \left( \frac{ P(\text{Robo}) }{ P(\text{Sin siniestro}) } \right), \]

manteniendo constantes las demás variables.

---

9.8 Exponenciación de los coeficientes

En el modelo multinomial también podemos calcular:

\[ e^{\beta}. \]

Estas cantidades expresan cambios multiplicativos en los odds relativos de una categoría respecto de la categoría base.

rrr_multi <- exp(
  coef(
    modelo_multinomial
  )
)

rrr_multi

##                  (Intercept)   edad siniestros_previos zonaMedia zonaAlta
## Danos materiales     0.12834 0.9881              1.681     1.359    1.480
## Lesiones             0.06550 0.9719              1.631     1.576    1.709
## Robo                 0.03886 0.9921              1.303     1.884    3.885
##                  uso_vehiculoComercial valor_asegurado
## Danos materiales                1.4075               1
## Lesiones                        2.2520               1
## Robo                            0.9271               1

Organizamos la información.

tabla_rrr_multi <- as.data.frame(
  rrr_multi
) %>%
  rownames_to_column(
    "Categoria"
  ) %>%
  pivot_longer(
    cols = -Categoria,
    names_to = "Variable",
    values_to = "RRR"
  )

knitr::kable(
  tabla_rrr_multi,
  digits = 4,
  caption = "Razones exponenciadas del modelo multinomial"
)

Razones exponenciadas del modelo multinomial

Categoria	Variable	RRR
Danos materiales	(Intercept)	0.1283
Danos materiales	edad	0.9881
Danos materiales	siniestros_previos	1.6805
Danos materiales	zonaMedia	1.3586
Danos materiales	zonaAlta	1.4801
Danos materiales	uso_vehiculoComercial	1.4075
Danos materiales	valor_asegurado	1.0000
Lesiones	(Intercept)	0.0655
Lesiones	edad	0.9719
Lesiones	siniestros_previos	1.6307
Lesiones	zonaMedia	1.5756
Lesiones	zonaAlta	1.7094
Lesiones	uso_vehiculoComercial	2.2520
Lesiones	valor_asegurado	1.0000
Robo	(Intercept)	0.0389
Robo	edad	0.9921
Robo	siniestros_previos	1.3029
Robo	zonaMedia	1.8838
Robo	zonaAlta	3.8854
Robo	uso_vehiculoComercial	0.9271
Robo	valor_asegurado	1.0000

---

9.9 Ejemplo: zona alta y robo

Extraemos:

rrr_robo_zona_alta <-
  rrr_multi[
    "Robo",
    "zonaAlta"
  ]

rrr_robo_zona_alta

## [1] 3.885

El valor es aproximadamente:

\[ 3.885. \]

9.9.1 Interpretación actuarial

Manteniendo constantes las demás variables, pertenecer a una zona alta multiplica la razón:

\[ \frac{ P(\text{Robo}) }{ P(\text{Sin siniestro}) } \]

por aproximadamente:

\[ 3.89 \]

respecto de una póliza localizada en zona baja.

Nuevamente:

\[ e^\beta \]

no debe interpretarse automáticamente como una multiplicación equivalente de la probabilidad.

---

9.10 Predicción de probabilidades multinomiales

Podemos obtener una probabilidad para cada categoría.

probabilidades_multi <- predict(
  
  modelo_multinomial,
  
  newdata =
    seguros_multi,
  
  type = "probs"
  
)

head(
  probabilidades_multi,
  10
)

##    Sin siniestro Danos materiales Lesiones    Robo
## 1         0.7917          0.08011  0.01351 0.11470
## 2         0.7466          0.13526  0.06311 0.05502
## 3         0.7658          0.14104  0.03115 0.06206
## 4         0.8333          0.08707  0.01775 0.06184
## 5         0.7569          0.12118  0.04055 0.08139
## 6         0.7284          0.14894  0.03057 0.09211
## 7         0.7619          0.13885  0.02474 0.07451
## 8         0.6767          0.14228  0.06894 0.11204
## 9         0.6695          0.21247  0.06718 0.05090
## 10        0.8076          0.12166  0.02677 0.04396

Cada fila contiene:

\[ P(\text{Sin siniestro}), \]

\[ P(\text{Daños materiales}), \]

\[ P(\text{Lesiones}), \]

\[ P(\text{Robo}). \]

Comprobamos que sumen 1.

head(
  rowSums(
    probabilidades_multi
  ),
  10
)

##  1  2  3  4  5  6  7  8  9 10 
##  1  1  1  1  1  1  1  1  1  1

---

9.11 Predicción para perfiles actuariales

Construiremos tres perfiles.

perfiles_multi <- tibble(
  
  perfil = c(
    "Perfil A",
    "Perfil B",
    "Perfil C"
  ),
  
  edad = c(
    55,
    38,
    24
  ),
  
  siniestros_previos = c(
    0,
    1,
    3
  ),
  
  zona = factor(
    c(
      "Baja",
      "Media",
      "Alta"
    ),
    levels = levels(
      seguros_multi$zona
    )
  ),
  
  uso_vehiculo = factor(
    c(
      "Particular",
      "Particular",
      "Comercial"
    ),
    levels = levels(
      seguros_multi$uso_vehiculo
    )
  ),
  
  valor_asegurado = c(
    90000,
    140000,
    220000
  )
)

prob_perfiles_multi <- predict(
  
  modelo_multinomial,
  
  newdata =
    perfiles_multi,
  
  type = "probs"
  
)

resultado_perfiles_multi <- bind_cols(
  perfiles_multi,
  as_tibble(
    prob_perfiles_multi
  )
)

knitr::kable(
  resultado_perfiles_multi,
  digits = 4,
  caption = "Probabilidades multinomiales para tres perfiles actuariales"
)

Probabilidades multinomiales para tres perfiles actuariales

perfil	edad	siniestros_previos	zona	uso_vehiculo	valor_asegurado	Sin siniestro	Danos materiales	Lesiones	Robo
Perfil A	55	0	Baja	Particular	90000	0.8884	0.0677	0.0134	0.0305
Perfil B	38	1	Media	Particular	140000	0.7083	0.1631	0.0472	0.0814
Perfil C	24	3	Alta	Comercial	220000	0.2785	0.3716	0.1964	0.1535

9.11.1 Interpretación

Para cada asegurado obtenemos un vector:

\[ [ P(\text{Sin siniestro}), P(\text{Daños}), P(\text{Lesiones}), P(\text{Robo}) ]. \]

Esto permite distinguir entre asegurados que podrían presentar probabilidades globales de siniestro similares, pero distinta composición del riesgo.

Por ejemplo:

• un asegurado puede concentrar su riesgo en daños materiales;
• otro puede concentrarlo en robo;
• otro puede presentar una mayor probabilidad relativa de lesiones.

Esta información puede ser relevante porque los diferentes tipos de siniestro presentan:

• severidades distintas;
• costos de reparación distintos;
• distintos deducibles;
• distintas necesidades de prevención;
• diferentes patrones territoriales.

---

9.12 Clasificación multinomial

También podemos identificar la categoría con mayor probabilidad.

seguros_multi <- seguros_multi %>%
  mutate(
    
    clase_predicha = predict(
      modelo_multinomial,
      newdata = seguros_multi,
      type = "class"
    )
    
  )

matriz_multi <- table(
  Observado =
    seguros_multi$tipo_siniestro,
  Predicho =
    seguros_multi$clase_predicha
)

matriz_multi

##                   Predicho
## Observado          Sin siniestro Danos materiales Lesiones Robo
##   Sin siniestro             5308                4        0    0
##   Danos materiales           895                6        0    0
##   Lesiones                   273                0        0    0
##   Robo                       511                2        0    1

9.12.1 Interpretación

La matriz compara:

• categoría real;
• categoría predicha.

Sin embargo, debe evitarse evaluar el modelo exclusivamente mediante exactitud.

En seguros normalmente existe una clase dominante:

\[ P(\text{Sin siniestro}) > P(\text{cada tipo particular de siniestro}). \]

Un modelo que predijera siempre la categoría mayoritaria podría obtener una exactitud aparentemente elevada y ser poco útil para detectar eventos importantes.

---

10 Caso actuarial integrado

Ahora utilizaremos el modelo binomial para analizar un asegurado específico.

Supongamos las siguientes características:

• edad: 30 años;
• antigüedad de la póliza: 3 años;
• siniestros previos: 2;
• zona: Alta;
• uso del vehículo: Particular;
• valor asegurado: 150,000.

nuevo_asegurado <- tibble(
  
  edad = 30,
  
  antiguedad_poliza = 3,
  
  siniestros_previos = 2,
  
  zona = factor(
    "Alta",
    levels = levels(
      entrenamiento$zona
    )
  ),
  
  uso_vehiculo = factor(
    "Particular",
    levels = levels(
      entrenamiento$uso_vehiculo
    )
  ),
  
  valor_asegurado = 150000
  
)

prob_nuevo <- predict(
  
  modelo_logit,
  
  newdata =
    nuevo_asegurado,
  
  type = "response"
  
)

prob_nuevo

##      1 
## 0.3525

10.0.1 Interpretación

La probabilidad estimada es:

\[ \hat p= 0.3525. \]

En porcentaje:

\[ 35.25%. \]

Esto significa que, dadas las características introducidas y la estructura del modelo, la probabilidad estimada de presentar al menos un siniestro durante el período es aproximadamente:

\[ 35.25%. \]

---

10.1 ¿Qué ocurre si eliminamos los siniestros previos?

Mantendremos exactamente las mismas características, pero cambiaremos:

\[ SiniestrosPrevios=2 \]

por:

\[ SiniestrosPrevios=0. \]

asegurado_sin_historial <-
  nuevo_asegurado %>%
  mutate(
    siniestros_previos = 0
  )

prob_sin_historial <- predict(
  
  modelo_logit,
  
  newdata =
    asegurado_sin_historial,
  
  type = "response"
  
)

comparacion_historial <- tibble(
  
  Escenario = c(
    "2 siniestros previos",
    "0 siniestros previos"
  ),
  
  Probabilidad = c(
    as.numeric(
      prob_nuevo
    ),
    as.numeric(
      prob_sin_historial
    )
  )
  
)

knitr::kable(
  comparacion_historial,
  digits = 4,
  caption = "Efecto del historial de siniestros sobre la probabilidad estimada"
)

Efecto del historial de siniestros sobre la probabilidad estimada

Escenario	Probabilidad
2 siniestros previos	0.3525
0 siniestros previos	0.1213

10.1.1 Interpretación

Esta comparación ilustra una ventaja fundamental de un modelo multivariable.

Podemos modificar una característica y mantener las demás constantes.

Esto permite estudiar la asociación del historial de siniestros sin confundirla con diferencias de:

• edad;
• zona;
• valor asegurado;
• antigüedad;
• uso.

---

10.2 ¿Qué ocurre si cambiamos únicamente la zona?

asegurado_zona_baja <-
  nuevo_asegurado %>%
  mutate(
    
    zona = factor(
      "Baja",
      levels = levels(
        entrenamiento$zona
      )
    )
    
  )

prob_zona_baja <- predict(
  
  modelo_logit,
  
  newdata =
    asegurado_zona_baja,
  
  type = "response"
  
)

comparacion_zona <- tibble(
  
  Zona = c(
    "Alta",
    "Baja"
  ),
  
  Probabilidad = c(
    as.numeric(
      prob_nuevo
    ),
    as.numeric(
      prob_zona_baja
    )
  )
  
)

knitr::kable(
  comparacion_zona,
  digits = 4,
  caption = "Comparación de la probabilidad estimada según zona"
)

Comparación de la probabilidad estimada según zona

Zona	Probabilidad
Alta	0.3525
Baja	0.1959

10.2.1 Interpretación

Aunque el Odds Ratio de una variable es constante dentro del modelo, el cambio absoluto en probabilidad depende del perfil inicial.

Por ello:

\[ \Delta p \]

no es constante para todos los asegurados.

Esta es una consecuencia directa de la no linealidad de la función logística.

---

11 Errores frecuentes de interpretación

11.1 Error 1: interpretar \(\beta\) como cambio en probabilidad

Incorrecto:

Un coeficiente de 0.50 significa que la probabilidad aumenta 50%.

Correcto:

\[ \beta=0.50 \]

significa que los log-odds aumentan 0.50.

El Odds Ratio es:

\[ e^{0.50}\approx1.65. \]

---

11.2 Error 2: interpretar el Odds Ratio como Risk Ratio

Incorrecto:

OR = 2 significa que la probabilidad se duplica.

Correcto:

OR = 2 significa que los odds se duplican.

---

11.3 Error 3: utilizar únicamente p-valores

Un modelo no debe evaluarse únicamente porque sus coeficientes sean estadísticamente significativos.

También debemos estudiar:

• magnitud de efectos;
• estabilidad;
• AUC;
• calibración;
• coherencia actuarial.

---

11.4 Error 4: interpretar el pseudo \(R^2\) como el \(R^2\) lineal

Incorrecto:

Un pseudo \(R^2=0.20\) significa que explicamos 20% de la varianza.

Esta interpretación no es válida.

---

11.5 Error 5: utilizar AUC como medida de calibración

El AUC mide:

\[ \text{discriminación}. \]

No mide:

\[ \text{calibración}. \]

---

11.6 Error 6: evaluar únicamente la exactitud

Con clases desbalanceadas, una elevada exactitud puede ser engañosa.

---

---

12 Síntesis conceptual

La regresión logística binomial modela:

\[ Y\in\{0,1\}. \]

Su ecuación fundamental es:

\[ \boxed{ \log \left( \frac{p}{1-p} \right) = X^\top\beta } \]

La probabilidad se recupera mediante:

\[ \boxed{ p= \frac{ e^{X^\top\beta} }{ 1+e^{X^\top\beta} } } \]

Los coeficientes se transforman en Odds Ratios utilizando:

\[ \boxed{ OR=e^\beta } \]

La evaluación de un modelo debe considerar varias dimensiones.

12.0.1 Inferencia

Pregunta:

¿Existe evidencia de asociación entre las variables y el riesgo?

Herramientas:

• deviance;
• Wald;
• p-valores;
• intervalos de confianza.

12.0.2 Magnitud

Pregunta:

¿Cuánto cambian los odds?

Herramienta:

\[ OR=e^\beta. \]

12.0.3 Ajuste

Pregunta:

¿El modelo mejora respecto de una especificación más simple?

Herramientas:

• log-verosimilitud;
• pseudo \(R^2\);
• AIC.

12.0.4 Discriminación

Pregunta:

¿El modelo separa correctamente riesgos altos y bajos?

Herramientas:

• curva ROC;
• AUC;
• sensibilidad;
• especificidad.

12.0.5 Calibración

Pregunta:

¿Las probabilidades pronosticadas corresponden con las frecuencias observadas?

Herramientas:

• gráfico de calibración;
• Hosmer-Lemeshow;
• razón observado/esperado.

---

13 Aplicación dentro de un marco actuarial

Podemos visualizar el proceso de la siguiente manera:

\[ \boxed{ \text{Características del riesgo} } \]

\[ \downarrow \]

\[ \boxed{ \eta=X^\top\beta } \]

\[ \downarrow \]

\[ \boxed{ p= \frac{e^\eta}{1+e^\eta} } \]

\[ \downarrow \]

\[ \boxed{ \text{Probabilidad estimada de ocurrencia} } \]

\[ \downarrow \]

\[ \boxed{ \text{Segmentación y decisiones actuariales} } \]

Una regresión logística puede aportar información para:

• tarificación;
• suscripción;
• gestión de cartera;
• prevención;
• retención;
• fraude;
• clasificación de riesgos.

Sin embargo, un modelo estadístico no sustituye:

• criterio actuarial;
• validación;
• análisis de estabilidad;
• restricciones regulatorias;
• revisión de calidad de datos;
• análisis económico de las decisiones.

---

14 Conclusiones

La regresión logística permite convertir una combinación de características observables en una probabilidad de ocurrencia.

En Ciencias Actuariales esta propiedad resulta especialmente importante porque gran parte de la actividad profesional consiste en cuantificar incertidumbre.

El análisis no debe limitarse a preguntar:

¿El coeficiente es significativo?

También debemos preguntar:

¿El efecto es actuarialmente relevante?

¿El modelo discrimina adecuadamente?

¿Las probabilidades están calibradas?

¿El modelo mantiene su desempeño fuera de la muestra utilizada para estimarlo?

¿La especificación es coherente con el fenómeno asegurador?

Por tanto, un modelo logístico profesional debería ser evaluado considerando:

\[ \boxed{ \text{Coherencia actuarial} + \text{Interpretabilidad} + \text{Inferencia} + \text{Discriminación} + \text{Calibración} + \text{Validación} } \]

La regresión multinomial extiende esta lógica a situaciones en las que existen múltiples categorías de resultado, permitiendo estudiar de forma diferenciada distintos tipos de siniestro.

El objetivo final no es únicamente producir una clasificación, sino generar conocimiento cuantitativo sobre la estructura del riesgo.

---

15 Código de consulta rápida

A continuación se presenta una síntesis de las principales instrucciones utilizadas durante la clase.

15.1 Ajustar regresión logística

modelo <- glm(
  y ~ x1 + x2 + x3,
  data = datos,
  family = binomial(
    link = "logit"
  )
)

15.2 Resumen

summary(
  modelo
)

15.3 Probabilidades predichas

predict(
  modelo,
  type = "response"
)

15.4 Odds Ratios

exp(
  coef(
    modelo
  )
)

15.5 Intervalos para Odds Ratios

beta <- coef(
  modelo
)

se <- sqrt(
  diag(
    vcov(
      modelo
    )
  )
)

exp(
  cbind(
    beta - 1.96 * se,
    beta + 1.96 * se
  )
)

15.6 Comparación mediante deviance

anova(
  modelo_nulo,
  modelo,
  test = "Chisq"
)

15.7 AUC

roc_obj <- roc(
  y,
  probabilidad
)

auc(
  roc_obj
)

15.8 Hosmer-Lemeshow

hoslem.test(
  y,
  probabilidad,
  g = 10
)

15.9 Regresión multinomial

modelo_multi <- multinom(
  y ~ x1 + x2 + x3,
  data = datos
)

15.10 Probabilidades multinomiales

predict(
  modelo_multi,
  type = "probs"
)