Definiciones

La regresión logística es un método de regresión que permite estimar la probabilidad de una variable dicotómica en función de una o mas covariables. Se utiliza principalmente para tareas de clasificación.

Conceptualmente se trata de un Modelo Lineal Generalizado donde se utilizan las covariables para modelar el parámetro de la distribución de la variable a predecir.

Cuando la variable respuesta es binaria (toma valores 0 o 1 para cada observación individual), se modela mediante una distribución de Bernoulli, que es un caso especial de la Binomial con \(n = 1\):

\(Y \sim {\sf Bernoulli}(p)\)

\(Y = \begin{cases} 1 \;\;\;\; p \\ 0 \;\;\;\; 1-p \end{cases}\)

En este marco, la distribución tiene un único parámetro libre a estimar: la probabilidad del evento \(p\).

Modelamos la probabilidad del evento (la \(p\) de \(Y\)) en función de las covariables \((x_1, x_2, x_3....x_n)\).

Como los valores que puede tomar la probabilidad del evento se ubican entre 0 y 1, debemos interponer una función que permita mapear los números reales (que son los resultados del componente lineal \(\beta_0\) + \(\beta_1\) \(X_1\) + …. \(\beta_p\) \(X_p\)) en el espacio entre 0 y 1. Esta función interpuesta se conoce como función \(Link\) y la que mas frecuentemente se utiliza es la función \(logit\).

\(logit(P(Y = 1|X)) =\) \(\beta_0\) + \(\beta_1\) \(X_1\) + …. \(\beta_p\) \(X_p\)

\(logit(t) =\) \(\log(\frac{t}{1-t})\)

\(\frac{P(Y = 1|X)}{1-P(Y = 1|X)} =\) \(e^{\beta_0 + \beta_1 X_1 + .... \beta_p X_p}\)

El \(logit\) es el \(logaritmo\; natural\) afectando al \(odds\), no a la probabilidad directamente, por lo que la relación lineal existe entre \(X\) y el \(log(odds)\) del evento.

Resumen:

  • Los modelos lineales generalizados, de los que forma parte la regresión logística, utilizan una función para establecer la relación lineal entre las covariables y los parámetros de la distribución de la variable a predecir.

  • En regresión lineal, se modela el valor promedio de \(Y\) en función de \(X\). En Logística, se modela la \(probabilidad\) de que \(Y\) pertenezca a la clase \(1\) utilizando \(odds\) como medida de probabilidad. A estos \(odds\) se les aplica una función \(link\) llamada \(logit\) para permitir que la relación entre \(X\) y la probabilidad de que \(Y\) pertenezca a la clase \(1\) sea \(LINEAL\).

  • Relación \(LINEAL\) entre \(log(odds)\) y \(X\).

Ajuste

Una vez obtenida la relación lineal entre el logaritmo de los \(odds\) y la variable predictora \(X\), se tienen que estimar los parámetros \(\beta_0\) y \(\beta_1\). La combinación óptima de valores será aquella que tenga la máxima verosimilitud (\(maximum\; likelihood\)), es decir el valor de los parámetros \(\beta_0\) y \(\beta_1\) con los que se maximiza la probabilidad de obtener los datos observados. En el caso de la regresión multiple es igual pero estimando \(\beta_0, \beta_1, \beta_2....\beta_n\).

Los coeficientes \(\beta\) son desconocidos, y deben ser estimados basados en los datos con los que se entrena el modelo. En el caso de la regresión lineal se estiman utilizando la técnica de mínimos cuadrados. En el caso de otros \(GLM\) se utiliza el método de \(Maximum \;Likelihood\).

Se buscan iterativamente valores \(\beta_0, \beta_1, \beta_2....\beta_n\), tales que la probabilidad predicha de \(Y\) para cada observación individual, corresponda lo máximo posible con la realidad observada \((0\;o\;1)\). Dicho de otra forma, se buscan los valores \(\beta_0, \beta_1, \beta_2....\beta_n\) tales que, al utilizarlos en la fórmula \(p(X) = \frac{e^{\beta_0 + \beta_1 X_1 + .... \beta_p X_p}}{1+e^{\beta_0 + \beta_1 X_1 + .... \beta_p X_p}}\), se obtenga un valor cercano a \(0\) para los que tienen valor observado \(Y = 0\) y cercano a \(1\) para los que tienen valor observado \(Y=1\).

La contribución de cada par \((x_i, y_i)\) a la función de verosimilitud puede expresarse de forma compacta como:

\(\pi(x_i)^{y_i} [1-\pi(x_i)]^{1-y_i}\)

De manera que si la predicción resultó en \(1\), el término que aportará a la fórmula será \(\pi(x_i)^{y_i}\), y si fuera \(0\), \([1-\pi(x_i)]^{1-y_i}\). Dado que las observaciones se asumen \(independientes\), la \(función\; de\; verosimilitud\) se obtiene del producto de los términos anteriores para todas las observaciones:

\(l(\beta) = \prod_{i = 1}^{n}\pi(x_i)^{y_i} [1-\pi(x_i)]^{1-y_i}\)

El principio de \(máxima \;verosimilitud\) indica que el coeficiente \(\beta\) seleccionado será aquel que maximice la ecuación anterior. Matemáticamente, resulta mas sencillo trabajar con el \(log-likelihood\) por lo que se buscará maximizar \(L(\beta) = \ln[l(\beta)]\)

\(L(\beta_0,\beta_1,\beta_2....\beta_n) = \sum_{i=1}^{n} Y_i(\beta_0+\beta_1 x_i)-\sum_{i=1}^{n} \ln(1+e^{(\beta_0+\beta_1 x_i)})\)

Entonces el \(estimador\;de\;máxima\;verosimilitud\) (el valor de \(\beta\) elegido) será aquel que maximice el valor del \(log-likelihood\) en la \(función\; de\; verosimilitud\). Esta función no se resuelve analíticamente. Se debe iterar hasta maximizar el valor de \(L(\beta)\). Las técnicas de iteración utilizadas pueden ser Gradient Descent, Newton Raphson y otras.

La estimación por mínimos cuadrados utilizada en regresión lineal es un caso especial de \(máxima\; verosimilitud\) (válido cuando se asume distribución Normal de los errores).

Evaluación

Los Coeficientes

Cuando nos proponemos testear la significancia de alguna variable en cualquier modelo, en realidad nos estamos haciendo la siguiente pregunta: \("El\; modelo \;que\; incluye\; a\; esta\; variable,\; nos\; dice\; mas\; sobre\; la\; variable\; outcome\) \(que\; un\; modelo\; que\; NO\; la\; incluye?"\). Esta pregunta puede responderse comparando los valores observados de la variable outcome con los predichos por cada uno de los modelos (el que tiene y el que no tiene la variable en cuestión). Es importante recordar que en esta instancia, no estamos evaluando si los valores predichos se ajustan a lo observado (eso se llama \(Bondad\; de\; Ajuste\) y sera estudiado mas adelante).

La evaluación de significancia estadística de cada variable regresora se ilustra facilmente en regresion lineal. Cuando se evalua la significancia de los coeficientes utilizando una Tabla de ANOVA se pueden observar:

  1. La suma de los cuadrados de las distancias entre los residuos y el modelo nulo representado por una linea de pendiente \(0\) ubicado a la altura del promedio de \(Y\) (\(\beta_0\))

\(SST = \sum_{i=1}^{n} (y_i - \overline{y})^2\)

  1. La suma de los cuadrados de las distancias entre los residuos y el modelo que incluye a la variable.

\(SSRes = \sum_{i=1}^{n} (y_i - \hat{y})^2\) [Tambien se le dice \(SSE\)]

  1. La suma de los cuadrados de las distancias entre el modelo ajustado y el modelo nulo.

\(SSM = \sum_{i=1}^{n} (\hat{y_i} - \overline{y})^2\)

\(SST = SSM + SSRes\)

En el modelo que no contiene la variable regresora que esta siendo evaluada, el único parámetro es \(\beta_0\) que corresponde a \(\overline{Y}\), el promedio de la variable outcome. Cuando incluimos la variable en cuestión al modelo, cualquier descenso en el valor de \(SSRes\) respecto de \(SST\) se deberá al hecho de que la pendiente/coeficiente (\(\beta\)) de la variable es diferente de \(0\). En este caso se estudia la diferencia existente entre \(SST\) y \(SSRes\); mientras mayor sea, mayor relevancia tendrá la variable para explicar el outcome.

Conceptualmente la evaluación que se realiza en el modelo logístico es igual: se comparan los valores reales de la variable outcome con las predicciones del modelo \(con\) y \(sin\) la variable en cuestión.

Para esta comparación se utiliza la \(Deviance\), considerando que los valores observados son equivalentes a los predichos por un modelo saturado.

\(D = -2\ln\left[\frac{l(\hat{\beta})}{l_{saturado}}\right]\)

Donde \(l(\hat{\beta})\) es la verosimilitud del modelo ajustado y \(l_{saturado}\) es la verosimilitud del modelo saturado (aquel que tiene tantos parámetros como observaciones y reproduce exactamente los datos). Desarrollando:

\(D = -2\ln[l(\hat{\beta})] - (-2\ln[l_{saturado}])\)

Como el modelo saturado ajusta los datos perfectamente, su log-verosimilitud se utiliza como cota de referencia. En la práctica, para datos binarios individuales el modelo saturado asigna probabilidad 1 a cada evento y 0 a cada no-evento, con lo que \(\ln[l_{saturado}] = 0\), y entonces:

\(D = -2\ln[l(\hat{\beta})]\)

La relación \(l(\hat{\beta})/l_{saturado}\) se llama \(likelihood\;ratio\) y se multiplica su logaritmo natural por \(-2\) para obtener una cantidad cuya distribución sea conocida (chi-cuadrado).

El estadístico \(D\) o \(Deviance\) juega el mismo rol en regresión logística que el \(SSRes\) lo hace en regresión lineal.

Para evaluar la significancia de una variable independiente se compara el valor de \(D\) con y sin la variable en cuestión. Esta prueba se llama \(likelihood\;ratio\;test\):

\(G=D(modelo\; sin\; la\; variable) - D(modelo\; con\; la\; variable)\)

\(G= -2\ln\left[\frac{l(modelo\;sin\;variable)}{l(modelo\;con\;la\;variable)}\right]\)

El estadístico \(G\) (\(likelihood\;ratio\;test\)) es equivalente al \(F-test\) de la regresión lineal.

Interpretación

A diferencia de la regresión lineal, en la que \(\beta_1\) se corresponde con el cambio promedio en la variable dependiente \(Y\) debido al incremento en una unidad del predictor \(X\), en regresión logística \(\beta_1\) indica el cambio en el logaritmo de \(ODDS\) debido al incremento de una unidad de \(X\), o lo que es lo mismo, multiplica los \(ODDS\) por \(e^{\beta_1}\). Dado que la relación entre \(p(Y)\) y \(X\) no es lineal, \(\beta_1\) no se corresponde con el cambio en la probabilidad de \(Y\) asociada con el incremento de una unidad de \(X\). Cuánto se incremente la probabilidad de \(Y\) por unidad de \(X\) depende del valor de \(X\), es decir, de la posición en la curva logística en la que se encuentre.

Covariables numéricas

Para una variable continua como la edad, el \(OR = e^{\hat{\beta}}\) se interpreta de la siguiente manera: por cada año adicional de edad, el \(odds\) del evento se multiplica por \(e^{\hat{\beta}}\). En la práctica, y aunque no sea matemáticamente riguroso, se suele expresar como: “por cada año que aumenta la edad, el odds del evento aumenta (o disminuye) un X%”, donde ese porcentaje corresponde a \((e^{\hat{\beta}} - 1) \times 100\). Esta aproximación es aceptable mientras el efecto sea moderado y la escala de la variable sea razonable.

Por ejemplo, si \(\hat{\beta}_{edad} = -0.03\), entonces \(OR = e^{-0.03} \approx 0.97\), lo que se leería como: “por cada año adicional de edad, el odds de supervivencia disminuye aproximadamente un 3%, manteniendo el resto de las variables constantes”.

Covariables categóricas

Para una variable categórica, cada nivel genera una variable \(dummy\) que se compara contra la categoría de referencia (el nivel basal). El \(OR = e^{\hat{\beta}}\) se interpreta como cuántas veces es mayor (o menor) el \(odds\) del evento en ese nivel respecto de la categoría de referencia, manteniendo el resto de las variables constantes.

Por ejemplo, si para la variable Pclass la primera clase es la referencia y \(\hat{\beta}_{Pclass3} = -1.8\), entonces \(OR = e^{-1.8} \approx 0.17\), lo que se leería como: “el odds de supervivencia en pasajeros de tercera clase es aproximadamente 0.17 veces el de los pasajeros de primera clase, es decir, un 83% menor”. Equivalentemente puede expresarse como: “el odds de no sobrevivir en tercera clase es aproximadamente 6 veces mayor que en primera clase”.

Predicción

Una vez que se conoce la relación entre los predictores y la variable outcome, se puede realizar la predicción utilizando una base de datos nueva o bien la partición destinada a testear el modelo construido, para estimar la validez externa del mismo.

Clasificación

Una de las principales aplicaciones de un modelo de regresión logística es clasificar la variable cualitativa en función del valor que tome el predictor. Para conseguir esta clasificación es necesario establecer un threshold de probabilidad a partir de la cual se considera que la variable pertenece a uno de los niveles \((0\;o\;1)\). Por ejemplo, se puede asignar una observación al grupo \(1\) si \(\hat{p}(Y = 1|X) > 0.5\) y al grupo \(0\) si no.

Evaluación del modelo

Performance del modelo evaluando los predictores

Se considera que el modelo es útil si es capaz de mostrar una mejora explicando las observaciones respecto al modelo nulo (sin predictores). El \(likelihood\;ratio\;test\) calcula la significancia de la diferencia de residuos entre el modelo de interés y el modelo nulo. El estadístico \(D\) sigue una distribución chi-cuadrado con grados de libertad equivalentes a la diferencia de grados de libertad de los dos modelos. Este test puede realizarse en \(R\) manualmente o bien generando una tabla de \(Anova\) para evaluar la \(Deviance\) del modelo agregando de a un predictor a la vez.

Calibración

La calibración evalúa el grado de concordancia entre las probabilidades predichas por el modelo y la frecuencia de eventos observados en los datos. Un modelo bien calibrado es aquel en el que, por ejemplo, entre los pacientes a quienes se les asigna una probabilidad predicha del 30%, aproximadamente el 30% efectivamente presenta el evento.

La evaluación de la calibración es independiente de la discriminación: un modelo puede discriminar bien (separar positivos de negativos) pero estar sistemáticamente sesgado en la escala de probabilidades, o viceversa.

Test de Hosmer-Lemeshow

El test de Hosmer-Lemeshow es una prueba de bondad de ajuste para modelos de clasificación binaria. Específicamente, el test calcula si la tasa de eventos observados (\(1s\)) coincide con la tasa de eventos predichos por el modelo. Para esto se separa a la población en cuantiles (tipicamente decilos) y se realiza una comparación basada en el estadístico \(\chi^2\). En este caso, dado que la hipótesis nula es “La tasa de eventos observados es igual a la tasa de eventos predicha”, si el \(p-valor\) del test resulta alto (zona de NO rechazo) estaría indicando que el modelo predice adecuadamente.

\(G_{HL}=\sum\limits_{j=1}^{g}\frac{(o_j -e_j)^2}{e_j}\)

Donde \(o_j\) son datos observados y \(e_j\) son datos predichos en cada grupo \(j\).

Esta prueba habitualmente genera una tabla comparando los decilos (o grupos seleccionados) y se reporta el estadístico \(G_{HL}\) con su \(p-valor\) asociado.

Limitaciones del test de Hosmer-Lemeshow: la prueba es sensible al tamaño muestral y al número de grupos elegido. Con muestras grandes, diferencias triviales en calibración pueden generar \(p-valores\) significativos; con muestras pequeñas, la prueba tiene escasa potencia. Adicionalmente, al resumir la calibración en un único estadístico global, puede no detectar problemas localizados en zonas específicas de la distribución de probabilidades.

Slope e Intercept de calibración

Un marco más informativo para evaluar la calibración es la regresión de calibración, que consiste en ajustar un modelo logístico donde el único predictor es el \(logit\) de las probabilidades predichas:

\(logit(Y) = \alpha + \beta \cdot logit(\hat{p})\)

En un modelo perfectamente calibrado se espera \(\alpha = 0\) y \(\beta = 1\).

  • El intercept de calibración (también denominado calibration-in-the-large) evalúa si la probabilidad media predicha coincide con la proporción de eventos observada. Un valor de \(\alpha \neq 0\) indica una sobreestimación o subestimación sistemática del riesgo.

  • El slope de calibración evalúa si el rango de predicciones es adecuado. Un valor \(\beta < 1\) es característico del sobreajuste (overfitting): el modelo genera predicciones extremas (muy cercanas a 0 o a 1) que no se corresponden con la realidad observada. Un valor \(\beta > 1\) indica que las predicciones son demasiado moderadas.

ICI y E50

Incluso cuando el slope y el intercept son adecuados, pueden existir problemas de calibración locales que esos estadísticos no capturan. El Integrated Calibration Index (ICI) y el E50 permiten una evaluación más granular.

Ambos se basan en estimar una curva de calibración suavizada (habitualmente mediante LOESS) que relaciona las probabilidades predichas con la proporción de eventos observados, y luego cuantificar cuánto se aleja esa curva de la diagonal perfecta (\(predicho = observado\)):

  • ICI: diferencia absoluta media entre la curva suavizada y la bisectriz, integrada sobre toda la distribución de probabilidades predichas. Puede interpretarse como el error de calibración promedio ponderado por la densidad de predicciones. Un ICI de 0 indica calibración perfecta.

  • E50: mediana de las diferencias absolutas entre la curva suavizada y la bisectriz. Es una medida robusta que refleja el error de calibración típico en el centro de la distribución de predicciones.

A diferencia del test de Hosmer-Lemeshow, estas métricas no dependen de la elección de grupos y evalúan la calibración en toda la escala de probabilidad, lo que las hace más sensibles y comparables entre estudios.

Performance del modelo para generar una discriminación adecuada

Una vez que uno tiene información sobre datos predichos y datos observados, se tiene suficiente material para generar una \(matriz\; de\; confusión\). Con ésta, podremos observar variaciones de la sensibilidad y especificidad del modelo corriendo el umbral o punto de corte a partir del cual consideraremos una predicción como \(1\;o\;0\). De esta forma podemos además generar una curva \(ROC\) para observar el comportamiento del modelo independientemente del punto de corte seleccionado.

El \(AUC\) (area bajo la curva \(ROC\)) hace referencia a la probabilidad de que, al seleccionar al azar un caso positivo (\(1\)) y uno negativo (\(0\)), el caso positivo tenga una probabilidad predicha mayor que el caso negativo. Un \(AUC = 0.5\) equivale a un modelo sin capacidad discriminativa (azar puro), mientras que \(AUC = 1\) corresponde a discriminación perfecta.

EJEMPLO

Definición de los datos

Se utilizará la base de datos de pasajeros del Titanic (train_titanic.csv), ampliamente conocida y disponible en Kaggle. El objetivo es predecir la supervivencia al hundimiento en función de características de los pasajeros.

Las variables utilizadas son:

Variable Definición Codificación
Survived Supervivencia 0 = No, 1 = Sí
Pclass Clase del pasaje 1 = Primera, 2 = Segunda, 3 = Tercera
Sex Sexo female / male
Age Edad en años continua
SibSp Número de hermanos/cónyuge a bordo continua
Fare Tarifa pagada continua 
library(tidyverse)
library(data.table)
library(sjPlot)
library(pROC)
library(vcdExtra)
library(ResourceSelection)
library(CalibrationCurves)
library(tableone)

options(scipen = 999)

# Carga y preparación
Titanic_full <- fread("train_titanic.csv") %>%
  dplyr::select(Survived, Pclass, Sex, Age, SibSp, Fare) %>%
  na.omit()

Titanic_full$Pclass   <- as.factor(Titanic_full$Pclass)
Titanic_full$Sex      <- as.factor(Titanic_full$Sex)
Titanic_full$Survived <- as.numeric(Titanic_full$Survived)

# Partición 70/30 para entrenamiento y validación
set.seed(123)
idx   <- sample(1:nrow(Titanic_full), 0.7 * nrow(Titanic_full))
train <- Titanic_full[idx, ]
test  <- Titanic_full[-idx, ]

Tabla descriptiva por outcome

catvars <- c("Pclass", "Sex")
vars    <- c("Pclass", "Sex", "Age", "SibSp", "Fare")

Tab_1 <- CreateTableOne(vars = vars, strata = "Survived",
                        factorVars = catvars, data = train)
kableone(Tab_1)
0 1 p test
n 265 171
Pclass (%) <0.001
1 45 (17.0) 65 (38.0)
2 49 (18.5) 47 (27.5)
3 171 (64.5) 59 (34.5)
Sex = male (%) 230 (86.8) 54 (31.6) <0.001
Age (mean (SD)) 30.62 (12.96) 27.87 (13.52) 0.034
SibSp (mean (SD)) 0.51 (1.21) 0.49 (0.74) 0.888
Fare (mean (SD)) 23.01 (31.03) 44.27 (60.05) <0.001

Construcción del modelo

El objetivo de esta sección no es la construcción paso a paso del modelo sino la comprensión de la interpretación de los coeficientes y de las métricas de performance.

Se construye un modelo con todas las variables disponibles para predecir la supervivencia.

mod <- glm(Survived ~ Pclass + Sex + Age + SibSp + Fare,
           family = "binomial", data = train)

t <- tab_model(mod, show.se = T, show.dev = T, show.intercept = F,
               show.loglik = T, digits = 2, digits.p = 3, show.r2 = F)

knitr::asis_output(t$knitr)
  Survived
Predictors Odds Ratios std. Error CI p
Pclass [2] 0.39 0.16 0.17 – 0.88 0.024
Pclass [3] 0.11 0.05 0.05 – 0.24 <0.001
Sex [male] 0.05 0.02 0.03 – 0.09 <0.001
Age 0.95 0.01 0.93 – 0.97 <0.001
SibSp 0.70 0.11 0.50 – 0.93 0.022
Fare 1.00 0.00 1.00 – 1.01 0.345
Observations 436
Deviance 377.238
log-Likelihood -188.619 
anova(mod, test = "LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Survived
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev              Pr(>Chi)    
## NULL                     435     584.00                          
## Pclass  2   40.202       433     543.80        0.000000001863 ***
## Sex     1  142.257       432     401.54 < 0.00000000000000022 ***
## Age     1   17.922       431     383.62        0.000023010228 ***
## SibSp   1    5.388       430     378.23               0.02028 *  
## Fare    1    0.992       429     377.24               0.31926    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Observamos el valor de los coeficientes y su \(p-valor\) asociado. En las variables categóricas podemos ver que se generó una variable \(dummy\) para cada nivel de la misma. Recordar que en este caso el coeficiente se refiere a variaciones en el \(log(odds)\) de supervivencia \(respecto\; de\; la\; clase\; basal\).

Respecto de las mediciones asociadas al modelo completo, observamos:

  • \(null.deviance:\) \(Deviance\) de un modelo nulo (sin predictores).
  • \(deviance:\) \(-2 \times logLik\) del modelo ajustado [se multiplica por \(-2\) para que sea mas intuitivo el hecho de buscar el valor mas bajo de \(deviance\)].
  • \(logLik:\) logaritmo de \(Likelihood\) del modelo.
  • \(AIC:\) Criterio de información de Akaike. Es una medida de evaluación del modelo que incluye evaluación de cantidad de parámetros utilizados y la \(Deviance\) del modelo: \(AIC = 2k - 2\ln(l(\hat{\beta}))\) donde \(k\) es la cantidad de parámetros del modelo. Recompensa la bondad de ajuste y penaliza el sobreajuste (\(overfitting\)). No se lo puede interpretar en valores absolutos sino que se utiliza para comparación entre modelos.

La tabla de \(ANOVA\) muestra la contribución a la disminución de la \(Deviance\) que tiene cada variable.

Interpretación de coeficientes

El coeficiente estimado para el \(intercept\) (\(\beta_0\)) es el valor esperado del \(logaritmo\; de\; odds\) de supervivencia con todas las variables continuas en \(cero\) y las categóricas en la categoría de referencia.

El coeficiente de Age indica que por cada año adicional de edad, el \(log(odds)\) de supervivencia disminuye en la magnitud de ese coeficiente, independientemente de las otras variables. En términos de \(ODDS\), esto equivale a multiplicarlos por \(e^{\hat{\beta}_{age}}\).

En el caso de la variable categórica Pclass, dado que involucra 3 niveles, los coeficientes se refieren a la modificación del \(log(odds)\) de supervivencia respecto del nivel de referencia (primera clase). Exponenciando los coeficientes se obtienen los \(Odds\;Ratios\) correspondientes.

Para Sex, el coeficiente asociado a male cuantifica la diferencia en el \(log(odds)\) de supervivencia respecto de las mujeres (categoría de referencia).

Calibración

La evaluación de calibración se realiza sobre el set de validación (test), es decir sobre observaciones que el modelo nunca vio durante el entrenamiento. Este es el escenario relevante para evaluar la utilidad clínica o predictiva del modelo.

# Predicciones en test
pred_test <- predict(mod, newdata = test, type = "response")

# --- Hosmer-Lemeshow ---
HL <- hoslem.test(test$Survived, pred_test, g = 10)

HL_DF <- data.frame(cbind(HL$observed, HL$expected)) %>%
  rownames_to_column(var = "quantilos") %>%
  dplyr::rename(observados = y1, esperados = yhat1) %>%
  dplyr::mutate(indice = as.factor(1:nrow(HL$observed)),
                CANT = y0 + observados) %>%
  dplyr::select(quantilos, CANT, observados, esperados, indice)

HL_DF %>% knitr::kable(caption = "Hosmer-Lemeshow: observados vs esperados por decil")
Hosmer-Lemeshow: observados vs esperados por decil
quantilos CANT observados esperados indice
[0.0205,0.0543] 19 1 0.7531914 1
(0.0543,0.0873] 19 5 1.3558579 2
(0.0873,0.0915] 19 1 1.7318552 3
(0.0915,0.125] 18 2 1.8984318 4
(0.125,0.244] 19 5 3.2065760 5
(0.244,0.413] 19 6 6.3562550 6
(0.413,0.582] 18 6 9.0558851 7
(0.582,0.724] 19 11 12.4279120 8
(0.724,0.883] 19 16 15.4751683 9
(0.883,0.986] 19 18 17.8828567 10 
HL
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  test$Survived, pred_test
## X-squared = 14.873, df = 8, p-value = 0.06166
ggplot(HL_DF, aes(x = indice)) +
  geom_point(aes(y = esperados, colour = "Esperados")) +
  geom_point(aes(y = observados, colour = "Observados")) +
  scale_colour_manual(values = c("Esperados" = "blue", "Observados" = "red")) +
  labs(x = "Decil de riesgo", y = "Proporción de eventos",
       colour = "", title = "Calibración por decilos - Set de validación")

## === Métricas de Calibración ===
## Intercept (Cal-in-the-large): 0.035
## Slope:                        0.826
## ICI:                          0.0389
## E50:                          0.045
  • Intercept (Cal-in-the-large): si es \(\approx 0\), las probabilidades medias predichas coinciden con la tasa de eventos observada.
  • Slope: si es \(\approx 1\), el rango de predicciones es adecuado. Un valor \(< 1\) indica que el modelo produce predicciones demasiado extremas (sobreajuste).
  • ICI: error de calibración promedio a lo largo de toda la escala de probabilidades predichas. Valores cercanos a 0 indican buena calibración.
  • E50: error de calibración mediano; refleja el error típico en el centro de la distribución de predicciones.

Un modelo bien calibrado mostrará los puntos de la curva alineados a lo largo de la diagonal, con ICI cercano a 0, slope \(\approx 1\) e intercept \(\approx 0\).

Clasificación - Discriminación

# Gráfico de violín: distribución de predichos según outcome real
df_pred <- data.frame(predichos = pred_test, observados = factor(test$Survived))

ggplot(df_pred, aes(observados, predichos)) +
  geom_violin(aes(colour = observados)) +
  geom_hline(yintercept = 0.5, linetype = "dashed") +
  labs(x = "Supervivencia observada", y = "Probabilidad predicha",
       title = "Distribución de probabilidades predichas por grupo")

# Curva ROC
roc_obj <- roc(test$Survived, pred_test)

plot(roc_obj, main = "Curva ROC - Set de validación",
     col = "steelblue", lwd = 2, legacy.axes = TRUE)
text(0.8, 0.2,
     paste("AUC =", round(auc(roc_obj), 3)),
     adj = c(0, 1))

En el gráfico de violín podemos evaluar cuán separadas están las distribuciones de probabilidades predichas entre quienes sobrevivieron y quienes no. Una buena discriminación se traduce en distribuciones con escasa superposición.

La curva \(ROC\) arroja un \(AUC\) que cuantifica la probabilidad de que, al seleccionar al azar un pasajero que sobrevivió y uno que no, el sobreviviente tenga asignada una probabilidad predicha mayor. Un \(AUC\) en el rango 0.70–0.80 suele considerarse aceptable; por encima de 0.80, bueno.

La evaluación de Calibración y Discriminación que vimos recién se realizó en el set de validación, que el modelo no utilizó para el entrenamiento. Esto permite una estimación más realista de la performance del modelo en datos nuevos.

Otros elementos utilizados en validación de modelos de clasificación como \(Cross-Validation\) o \(Bootstrapping\) no están incluidos en este documento.

Bibliografía

  1. James G, Witten D, Hastie T, Tibshirani R. An Introduction to Statistical Learning: with Applications in R, 2nd ed. Springer; 2021.

  2. Harrell FE. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis, 2nd ed. Springer; 2015.

  3. Austin PC, Steyerberg EW. The Integrated Calibration Index (ICI) and related metrics for quantifying the calibration of logistic regression models. Statistics in Medicine. 2019;38(21):4051–4065.

  4. Collins, G. S., Dhiman, P., Ma, J., Schlussel, M. M., Archer, L., Van Calster, B., … & Riley, R. D. (2024). Evaluation of clinical prediction models (part 1): from development to external validation. bmj, 384.

  5. Riley, R. D., Archer, L., Snell, K. I., Ensor, J., Dhiman, P., Martin, G. P., … & Collins, G. S. (2024). Evaluation of clinical prediction models (part 2): how to undertake an external validation study. bmj, 384.