3.6 Regresión logística

¡Atención! Aunque se llama regresión logística, la regresión logística es en realidad un método de clasificación.

Por Ariane Ducellier (2021)

🖥️ Diapositivas — Sesión 16 (mié 4 nov)

En este laboratorio vamos a hablar de:

• Un método de clasificación simple: la regresión logística
• El método de descenso de gradiente
• La diferenciación automática
• Una introducción a PyTorch

Regresión logística¶

Recuerde la regresión lineal:

$y = b + x w + \epsilon$

$y$ es un vector de longitud $n$, $x$ es una matriz con $n$ filas y $p$ columnas, correspondientes a $n$ observaciones y $p$ características que se usan para explicar $y$.

$b$ es un escalar. $w$ es un vector de longitud $p$. ε es un vector de error aleatorio, de longitud $n$. Es independiente de $x$ y tiene media cero.

Nuestro objetivo es encontrar los mejores valores de $b$ y $w$ de modo que los valores de $\hat{y} = b + x w$ queden tan cerca como sea posible de los valores reales $y$.

En la regresión lineal, $y$ es una variable cuantitativa. ¿Qué pasa si $y$ es una variable cualitativa, por ejemplo $y = 0$ para «no» y $y = 1$ para «sí»?

Una manera de usar la regresión para resolver un problema de clasificación es modelar la probabilidad de que la variable $y$ tome el valor 1:

$P (y = 1) = b + x w$

Una vez que hemos encontrado los mejores valores $\hat{b}$ y $\hat{w}$, calculamos $\hat{y} = \hat{b} + x \hat{w}$. Si $\hat{y} \geq 0.5$, decidimos clasificar esta observación como $y = 1$, es decir «sí». Si $\hat{y} < 0.5$, decidimos clasificar esta observación como $y = 0$, es decir «no».

Hay un problema con este método. Quisiéramos tener $0 \leq P (y = 1) \leq 1$ porque es una probabilidad. Sin embargo, nada en esta formulación obliga a $b$ y $w$ a tomar valores tales que $\hat{y} = \hat{b} + x \hat{w}$ caiga siempre en $[0, 1]$.

Para resolver este problema, podemos escribir en su lugar:

$z = b + x w$ y $P (y = 1) = \frac{1}{1 + e^{-z}}$

De esa manera, siempre tenemos $0 \leq P (y = 1) \leq 1$. Cuando $b + x w$ se hace grande, $P (y = 1)$ se acerca a 1, y el valor «sí» es cada vez más probable. Cuando $b + x w$ se hace pequeño, $P (y = 1)$ se acerca a 0, y el valor «no» es cada vez más probable.

¿Cómo encontramos los valores óptimos de $b$ y $w$? Definimos la función de pérdida de entropía cruzada. Para una observación, la función de pérdida es:

$\mathcal{L} = - \left(y_i \log \hat{y}_i + (1 - y_i) \log (1 - \hat{y}_i)\right)$ con $\hat{y}_i = \frac{1}{1 + e^{- (b + x_i^T w)}}$ donde $x_i$ es la $i$-ésima fila de x.

Si la observación verdadera $y_i$ es 1 («sí») y $\hat{y}_i = 1$, la función de pérdida vale 0. Si $\hat{y}_i = 0$, la función de pérdida tiende a infinito.

Si la observación verdadera $y$ es 0 («no») y $\hat{y}_i = 0$, la función de pérdida vale 0. Si $\hat{y}_i = 1$, la función de pérdida tiende a infinito.

Para las $n$ observaciones, escribimos:

$\mathcal{L} = \sum_{i = 1}^n \mathcal{L}_i$

Nuestro objetivo es entonces encontrar los valores de $b$ y $w$ que minimizan la función de pérdida. Note que, con esta formulación, $\mathcal{L}$ es siempre positiva.

Descenso de gradiente¶

Sabemos que el gradiente $\frac{\partial \mathcal{L}}{\partial w_j}$ es positivo si la pérdida $\mathcal{L}$ aumenta cuando $w_j$ aumenta. A la inversa, el gradiente $\frac{\partial \mathcal{L}}{\partial w_j}$ es negativo si la pérdida $\mathcal{L}$ disminuye cuando $w_j$ aumenta.

Para obtener valores cada vez más pequeños de la pérdida, en cada iteración tomamos:

$w_j^{(k + 1)} = w_j^{(k)} - \alpha \frac{\partial \mathcal{L}}{\partial w_j}$ para $j = 1 , \cdots , p$

$b^{(k + 1)} = b^{(k)} - \alpha \frac{\partial \mathcal{L}}{\partial b}$

Suponemos que el valor de α no es demasiado grande. Si el gradiente es positivo, el valor de $w_j$ disminuirá en cada iteración, y el valor de la función de pérdida disminuirá. Si el gradiente es negativo, el valor de $w_j$ aumentará en cada iteración, y el valor de la pérdida disminuirá.

Así que ahora todo lo que necesitamos hacer es calcular el gradiente de la función de pérdida.

Diferenciación automática¶

Hay tres maneras de calcular el gradiente. El primer método es usar la fórmula de la pérdida:

$\mathcal{L} (w_j , b) = - \sum_{i = 1}^n y_i \log (\frac{1}{1 + \exp (- b - \sum_{j = 1}^p w_j x_{i,j})}) + (1 - y_i) \log (1 - \frac{1}{1 + \exp (- b - \sum_{j = 1}^p w_j x_{i,j})})$

y calcular la fórmula exacta de las derivadas $\frac{\partial \mathcal{L}}{\partial w_j}$ y $\frac{\partial \mathcal{L}}{\partial b}$. Después solo hay que implementar la fórmula exacta en el código para calcular el gradiente.

Cuando la fórmula se complica más y más, se vuelve más y más probable cometer un error, ya sea en el cálculo de la fórmula de la derivada, ya sea en la implementación en el código.

El segundo método es calcular una aproximación del gradiente:

$\frac{\partial \mathcal{L}}{\partial w_j} = \frac{\mathcal{L}(w_j + \Delta w_j) - \mathcal{L}(w_j)}{\Delta w_j}$

Si se acumulan demasiadas aproximaciones, el método puede no funcionar muy bien y dar resultados inexactos.

El tercer método es usar la diferenciación automática. Si escribimos:

$z = x_i^T w + b = f_x(w, b)$, $\sigma = \frac{1}{1 + e^{-z}} = g(z)$ y $L = - (y_i \log(\sigma) + (1 - y_i) \log(1 - \sigma)) = h_y(\sigma)$, obtenemos:

$\frac{\partial L}{\partial w_j} = \frac{\partial f}{\partial w_j} g'(z) h'(\sigma)$

Es muy fácil calcular la fórmula exacta de las derivadas:

$\frac{\partial f}{\partial w_j}(w, b) = x_{i,j}$

$g'(z) = \frac{e^{-z}}{(1 + e^{-z})^2}$

$h'(\sigma) = - \frac{y_i}{\sigma} + \frac{1 - y_i}{1 - \sigma}$

Al calcular $L$, necesitamos entonces mantener en memoria los valores de $\frac{\partial f}{\partial w_j}(w, b)$, $g'(z)$ y $h'(\sigma)$ para poder calcular el gradiente. Eso es lo que hace PyTorch.

Introducción a PyTorch¶

PyTorch (https://pytorch.org/) es un paquete de Python que le permite construir y entrenar redes neuronales. Se basa en la diferenciación automática.

Los datos: potabilidad del agua¶

Usamos una tabla de mediciones de calidad del agua con un objetivo binario, Potability (0 = no potable, 1 = potable). Cada fila tiene 9 características numéricas (pH, dureza, sólidos disueltos, etcétera).

Algunas filas tienen valores faltantes. Las eliminamos y reiniciamos el índice de filas.

Divida antes de escalar¶

Apartamos el 30 % de las filas como conjunto de prueba antes de hacer cualquier otra cosa. Ajustamos el escalador solo sobre las filas de entrenamiento y lo aplicamos a ambos conjuntos: si el escalador viera las filas de prueba, información del conjunto de prueba se fugaría hacia el entrenamiento, y la evaluación final ya no sería honesta.

Antes de entrenar nada, calcule el modelo de referencia (baseline) trivial: un «modelo» que siempre predice la clase más común. Todo clasificador que valga la pena conservar debe superar este número.

Una observación, a mano¶

Vamos a calcular la pérdida correspondiente a la primera observación del conjunto de entrenamiento. En lugar de usar arreglos de Numpy para colocar nuestros datos y parámetros, vamos a usar tensores de torch, porque tienen propiedades que los arreglos de Numpy no tienen.

Estas son las características de la primera observación de entrenamiento:

Esta es la clase de la primera observación de entrenamiento:

Tomemos valores aleatorios para $w$ y $b$. Al crear estas variables usamos la opción requires_grad=True porque más adelante querremos calcular el gradiente con respecto a ellas.

Definamos $z = f(w, b) = x_i^T w + b$. Tenemos $\frac{\partial f}{\partial w_j} = x_{i,j}$ y $\frac{\partial f}{\partial b} = 1$.

Por defecto, PyTorch solo conserva los gradientes de los tensores hoja como W y B. Los tensores intermedios como $z$ no conservan sus gradientes; llamar a retain_grad() le pide a PyTorch almacenarlos para que podamos inspeccionarlos después de la pasada hacia atrás (backward pass).

Definamos $\sigma = g(z) = \frac{1}{1 + e^{-z}} = g(f(w, b)) = (g \circ f) (w, b)$. Tenemos $g'(z) = \frac{e^{-z}}{(1 + e^{-z})^2}$.

Note que aquí usamos la función torch.exp en lugar de numpy.exp. Eso es porque numpy solo calcula el valor de $e^x$ pero no sabe que la derivada de $e^x$ es $e^x$. Si queremos poder usar la diferenciación automática, necesitamos usar la función equivalente de torch, que calculará tanto $\sigma(z)$ como $\frac{\partial \sigma}{\partial z}(z)$. Este último valor será necesario cuando más adelante calculemos el gradiente.

Definamos $L = h(\sigma) = - (y_i \log(\sigma) + (1 - y_i) \log(1 - \sigma)) = h(g(z)) = h(g(f(w, b))) = (h \circ g \circ f) (w, b)$. Tenemos $L'(\sigma) = - (\frac{y_i}{\sigma} - \frac{1 - y_i}{1 - \sigma})$.

Ahora podemos calcular el gradiente de la pérdida para una observación. Este comando calcula el gradiente de L con respecto a todas las variables que conservan un gradiente, es decir W, B, z y sigma, pero no devuelve el valor.

Tenemos $\frac{\partial L}{\partial \sigma} = - (\frac{y_i}{\sigma} - \frac{1 - y_i}{1 - \sigma})$. Comparemos el resultado de PyTorch con la fórmula matemática exacta.

Tenemos $\frac{\partial L}{\partial z} = g'(z) h'(g(z))$, es decir $\frac{\partial L}{\partial z} = g'(z) h'(\sigma)$. Comparemos el resultado de PyTorch con la fórmula matemática exacta.

Tenemos $\frac{\partial L}{\partial b} = \frac{\partial f}{\partial b} g'(f(w, b)) h'(g(f(w, b)))$, es decir $\frac{\partial L}{\partial b} = \frac{\partial f}{\partial b} g'(z) h'(\sigma)$. Del mismo modo, tenemos $\frac{\partial L}{\partial w_j} = \frac{\partial f}{\partial w_j} g'(f(w, b)) h'(g(f(w, b)))$, es decir $\frac{\partial L}{\partial w_j} = \frac{\partial f}{\partial w_j} g'(z) h'(\sigma)$. Comparemos los resultados de PyTorch con las fórmulas matemáticas exactas.

Implementación de la regresión logística¶

Implementemos ahora la regresión logística usando todo el conjunto de entrenamiento.

Podríamos programar a mano la regla de actualización $w^{(k+1)} = w^{(k)} - \alpha \, \partial \mathcal{L} / \partial w$, y las ediciones anteriores de este libro lo hacían, con contabilidad de retain_grad en cada paso. torch.optim.SGD aplica exactamente esa regla de actualización, así que dejamos que se encargue de las actualizaciones de los parámetros y mantenemos explícita la fórmula de la pérdida. Escribimos la entropía cruzada binaria completa en lugar de llamar a torch.nn.BCELoss — el punto aquí es la transparencia, no la conveniencia.

Dos detalles prácticos:

• Acotamos σ lejos de 0 y de 1 antes de tomar logaritmos, para que la pérdida nunca devuelva log(0).
• Nos detenemos cuando el cambio relativo de la pérdida cae por debajo de 10^-6, o después de 2000 iteraciones.

Grafique la curva de pérdida. Debería decrecer rápido al principio y luego aplanarse.

Evaluación sobre el conjunto de prueba¶

Ahora predecimos tanto sobre el conjunto de entrenamiento como sobre el conjunto de prueba. Las ediciones anteriores de esta lección evaluaban el modelo sobre los mismos datos con los que fue entrenado; la brecha entre los dos números de abajo es exactamente la razón por la que eso exagera la habilidad.

Calculemos ahora algunas métricas de clasificación a mano, sobre el conjunto de prueba. N_test es el número de observaciones de prueba.

Ahora compare la exactitud de entrenamiento contra la exactitud de prueba, junto al modelo de referencia de clase mayoritaria.

La regresión logística apenas supera al modelo de referencia de clase mayoritaria en esta tabla, y la exhaustividad (recall) sobre la clase potable es baja. Ese es un resultado honesto: estas 9 características, combinadas linealmente, llevan poca señal sobre el objetivo. Un resultado casi nulo reportado contra un modelo de referencia es más útil que un score de entrenamiento inflado. Si alguien reporta solo una exactitud de entrenamiento, pregunte cuál era el modelo de referencia y qué dice un conjunto de prueba apartado.

¿Son honestas las probabilidades?¶

Todo lo anterior corta la salida de la sigmoide en el umbral 0.5 y califica las etiquetas resultantes. Pero el modelo produce una probabilidad, y en muchas aplicaciones la probabilidad misma es lo que se entrega. Un clasificador está calibrado cuando sus probabilidades declaradas coinciden con las frecuencias observadas: entre todas las muestras de prueba donde el modelo dice «70 % de probabilidad de potable», alrededor del 70 % deberían ser potables.

Dos herramientas lo miden sobre el conjunto de prueba:

• Un diagrama de confiabilidad (sklearn.calibration.calibration_curve) agrupa las muestras en intervalos según la probabilidad predicha y grafica la fracción observada de positivos en cada intervalo contra la probabilidad predicha media. Un modelo calibrado sigue la diagonal.
• El score de Brier es el error cuadrático medio de la probabilidad, $\frac{1}{N} \sum_{i=1}^N (\hat{y}_i - y_i)^2$. Más bajo es mejor. A diferencia de la exactitud, que solo ve de qué lado de 0.5 cae una predicción, el score de Brier castiga un 0.99 equivocado mucho más fuerte que un 0.55 equivocado.

Las probabilidades del modelo abarcan solo de 0.27 a 0.57 aproximadamente: nunca afirma certeza sobre ninguna muestra, lo cual es consistente con la señal débil de estas características. La curva oscila alrededor de la diagonal, con la mayoría de los intervalos centrales por debajo de ella — donde el modelo dice 40 % de potable, la frecuencia observada está más cerca de 30–40 %, una sobreconfianza leve hacia la clase potable. Cada intervalo por cuantiles contiene unas 60 muestras de prueba, así que cada fracción observada lleva un ruido de muestreo de aproximadamente ±0.06, y solo las desviaciones mayores que eso merecen interpretación. Veredicto: aproximadamente honesto, débilmente informativo. Las probabilidades son utilizables, pero en su mayor parte reformulan la tasa base.

Un competidor confiado, y cómo repararlo¶

Un método predict_proba no es garantía de que los números que salen de él merezcan llamarse probabilidades. Tome un bosque aleatorio pequeño: 10 árboles profundos, cada uno de los cuales casi memoriza el conjunto de entrenamiento. Su «probabilidad» de potable es la fracción de árboles que votan potable, y con 10 árboles sobreajustados esas fracciones caen con facilidad en valores extremos como 0.9 o 1.0.

El bosque memoriza el conjunto de entrenamiento y aun así supera al modelo logístico en exactitud de prueba. Ahora verifique si sus probabilidades merecen confianza, y repárelas si no. Lo que está en juego es práctico: cuando un modelo le dice a una oficina de gestión del riesgo «80 % de probabilidad», alrededor del 80 % de esos casos deberían cumplirse, o el número no es un pronóstico sobre el que alguien pueda actuar. La reparación reetiqueta los scores del modelo usando frecuencias observadas: aparte algunas muestras de datos, registre con qué frecuencia las muestras con score cercano a 0.9 resultan de verdad potables, y corrija cada score a esa frecuencia observada. Solo se permiten correcciones que preserven el orden — si el bosque puntuó la muestra A por encima de la muestra B, la probabilidad corregida de A queda igual o por encima de la de B —, de modo que el ordenamiento de las muestras sobrevive; solo cambian los números adjuntos. CalibratedClassifierCV proporciona las muestras apartadas reajustando el bosque sobre pliegues de validación cruzada (las divisiones rotativas de la lección 3.8), y el nombre que scikit-learn le da al ajuste que preserva el orden es regresión isotónica (method='isotonic').

La curva del bosque crudo cae muy por debajo de la diagonal a la derecha: donde afirma un 80 % o 90 % de probabilidad de potable, la frecuencia observada está cerca del 50 %. Miente con confianza exactamente en el rango sobre el que un usuario actuaría, y su exactitud de prueba no da ninguna pista de ello. La calibración isotónica baja el score de Brier y acerca la curva a la diagonal, aunque no la coloca sobre ella: el mapa de calibración se estima a partir de unas 1400 filas de entrenamiento y lleva su propio ruido de muestreo. La exactitud también mejora un poco, en parte porque CalibratedClassifierCV promedia los cinco bosques que ajusta a través de los pliegues — un pequeño bono de ensamble, no el punto del ejercicio.

En contextos de riesgo, la probabilidad es el entregable, no la etiqueta. Una ingeniera geotécnica reporta una probabilidad de licuefacción que entra en una verificación del reglamento de construcción; un gestor del agua actúa sobre una probabilidad de floración algal nociva; un mapa de llanura de inundación codifica probabilidades de excedencia de crecidas. La persona que fija el umbral de decisión — cuánta probabilidad justifica el costo de actuar — no suele ser la que entrenó el modelo, y tomará el número al pie de la letra. Un 90 % mal calibrado que en realidad significa 55 % traslada silenciosamente su error de modelado a la decisión de otra persona. Siempre que las probabilidades de un modelo alimenten un umbral que usted no controla, reporte un diagrama de confiabilidad y un score de Brier junto a la exactitud. El desacuerdo entre los miembros de un ensamble da una vista complementaria de la incertidumbre en la lección 3.9, y la misma verificación de confiabilidad se aplica a los ensambles profundos en la lección 4.5.

Apéndice¶

La regresión logística es un buen ejemplo para empezar a aprender sobre la diferenciación automática y PyTorch. Sin embargo, si de verdad quiere usar la regresión logística con su propio conjunto de datos, es mucho más fácil usar la función que ya existe en scikit-learn. Seguimos el mismo protocolo: ajustar sobre el conjunto de entrenamiento, reportar sobre el conjunto de prueba.

Leer los coeficientes¶

Las características fueron estandarizadas, así que cada coeficiente es el cambio en el log-odds (el logaritmo de la razón de probabilidades) de la potabilidad por un aumento de una desviación estándar en esa característica, manteniendo fijas las demás. Exponenciarlo lo convierte en una razón de probabilidades (odds ratio): un coeficiente $w_j$ multiplica la razón de probabilidades de potabilidad por $e^{w_j}$ por cada desviación estándar de la característica $j$. La celda de abajo lee de esta manera el coeficiente más grande.

Dos precauciones. Primera, el signo de un coeficiente describe el modelo ajustado, no la química del agua: con características correlacionadas y una señal tan débil, el ajuste puede asignarle un signo químicamente inverosímil a una característica mientras otra correlacionada absorbe el efecto contrario — la misma salvedad que la importancia de las características en la lección 3.7. Segunda, la LogisticRegression de sklearn aplica regularización L2 por defecto (C=1.0), así que sus coeficientes están encogidos hacia cero y difieren del ajuste sin penalización que entrenamos arriba por descenso de gradiente — solo ligeramente aquí, porque con unas 2000 muestras de entrenamiento la penalización por defecto es suave, pero la tabla de abajo muestra los dos vectores de coeficientes lado a lado para que usted verifique en lugar de suponer.