3.5 - Notas de Luis - Regresión logística desde cero¶
Cuatro meses después del proyecto de predicción de suministros, el Hospital Minermont se enfrenta a una nueva pregunta: ¿acudirá un paciente a su cita? En el conjunto del hospital se pierde alrededor del 18% de las citas, y en algunas especialidades hasta el 27%. Fisioterapia tiene la tasa más alta, así que el piloto empieza allí.
Esta vez no predecimos una cantidad, como las mascarillas del capítulo 2, sino una respuesta de sí o no: un problema de clasificación binaria. La regresión logística conserva la puntuación lineal de la regresión lineal, $z = w \cdot x + b$, y la pasa por la función sigmoide, que convierte cualquier número en una probabilidad entre 0 y 1:
$$\hat{p} = \sigma(z) = \frac{1}{1 + e^{-z}}$$
Si $\hat{p}$ supera un umbral (normalmente 0,5), predecimos que el paciente no acudirá.
En este cuaderno la construimos desde cero con NumPy, igual que hicimos con la regresión lineal:
- Importar las bibliotecas.
- Crear un conjunto de datos sintético de citas de Fisioterapia.
- Preparar los datos: conjuntos de entrenamiento y prueba, y escalado de características.
- Definir la función sigmoide.
- Definir la función de coste (entropía cruzada, o log loss).
- Entrenar el modelo con descenso del gradiente.
- Evaluarlo: matriz de confusión, sensibilidad, especificidad, VPP y VPN.
- Elegir el umbral de decisión.
- Interpretar los coeficientes.
- Predecir para pacientes nuevos.
Paso 1: Importar las bibliotecas¶
numpy para los cálculos, pandas para manejar los datos como una tabla y matplotlib para las gráficas. De scikit-learn solo usaremos una utilidad: dividir los datos en entrenamiento y prueba. El modelo lo escribimos nosotros.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split # solo para dividir los datos
np.set_printoptions(precision=3, suppress=True)
print("Bibliotecas cargadas.")
Bibliotecas cargadas.
Paso 2: Un conjunto de datos sintético de citas¶
Mientras Marta y Víctor reúnen los registros reales de asistencia del hospital, yo estoy creando un conjunto de datos sintético, como hicimos en el capítulo 2. Cada fila es una cita de Fisioterapia, descrita por las variables que el equipo consideró relevantes:
| Columna | Significado |
|---|---|
edad |
Edad del paciente (18–85 años) |
dias_desde_ultima_cita |
Días desde la cita anterior del paciente |
distancia_km |
Distancia de casa al hospital, en km (Minermont atiende una zona rural) |
hora |
Hora de la cita (de 8 a 19 h) |
primera_visita |
1 si es una primera visita, 0 si es una revisión |
ausencias_previas |
Número de citas a las que el paciente faltó anteriormente |
no_asiste |
Objetivo: 1 si el paciente no acudió, 0 si acudió |
La clase positiva (1) es la ausencia, porque es el suceso que queremos detectar. Detrás de los datos hay una regla oculta que solo usamos para generar las citas; el modelo tendrá que descubrir el patrón a partir de los ejemplos.
def generate_appointments(n, seed):
"""Citas sintéticas de Fisioterapia en Minermont: una fila por cita."""
rng = np.random.default_rng(seed)
age = rng.integers(18, 86, n) # años
days = rng.integers(1, 181, n) # días desde la cita anterior
distance = np.round(rng.gamma(2.0, 6.0, n).clip(0.5, 60), 1) # km hasta el hospital (zona rural)
hour = rng.integers(8, 20, n) # hora de la cita, de 8:00 a 19:00
first = rng.binomial(1, 0.35, n) # 1 = primera visita, 0 = revisión
previous = rng.poisson(0.6, n).clip(0, 6) # ausencias en el historial del paciente
late = np.where(hour >= 17, hour - 16, 0) # a última hora de la tarde se falta más
# La regla oculta detrás de los datos. Solo la usamos para generar las citas:
# el modelo nunca la ve y tiene que descubrir el patrón a partir de los ejemplos.
z = (-1.55 - 0.069 * (age - 50) + 0.026 * (days - 90) + 0.11 * (distance - 12)
+ 1.5 * (previous - 0.6) + 0.56 * (late - 0.75) + 0.63 * (first - 0.35))
no_show = rng.binomial(1, 1 / (1 + np.exp(-z))) # 1 = no acudió
return pd.DataFrame({
'edad': age, 'dias_desde_ultima_cita': days, 'distancia_km': distance, 'hora': hour,
'primera_visita': first, 'ausencias_previas': previous, 'no_asiste': no_show,
})
df = generate_appointments(1200, seed=42)
print(f"Citas: {len(df)}")
print(f"Tasa de ausencias: {df['no_asiste'].mean():.1%}")
df.head()
Citas: 1200 Tasa de ausencias: 26.2%
| edad | dias_desde_ultima_cita | distancia_km | hora | primera_visita | ausencias_previas | no_asiste | |
|---|---|---|---|---|---|---|---|
| 0 | 24 | 159 | 4.1 | 11 | 0 | 2 | 1 |
| 1 | 70 | 86 | 43.7 | 11 | 1 | 0 | 0 |
| 2 | 62 | 178 | 7.1 | 18 | 0 | 0 | 0 |
| 3 | 47 | 138 | 3.5 | 8 | 0 | 0 | 0 |
| 4 | 47 | 118 | 9.1 | 8 | 0 | 3 | 1 |
Antes de modelar, un vistazo rápido a los datos. ¿Cómo cambia la tasa de ausencias con algunas de las variables?
fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)
groups = [
(pd.cut(df['edad'], [17, 35, 50, 65, 85], labels=['18-35', '36-50', '51-65', '66-85']), 'Edad'),
(pd.cut(df['distancia_km'], [0, 5, 15, 30, 60], labels=['< 5', '5-15', '15-30', '> 30']), 'Distancia al hospital (km)'),
(df['ausencias_previas'].clip(upper=3).map({0: '0', 1: '1', 2: '2', 3: '3+'}), 'Ausencias previas'),
]
for ax, (group, title) in zip(axes, groups):
df.groupby(group, observed=True)['no_asiste'].mean().plot.bar(ax=ax, color='tab:blue')
ax.set_title(title)
ax.set_xlabel('')
ax.tick_params(axis='x', rotation=0)
axes[0].set_ylabel('Tasa de ausencias')
plt.tight_layout()
plt.show()
Los pacientes más jóvenes, los que viven lejos y los que ya han faltado a citas anteriores se ausentan más a menudo. Ninguna variable decide por sí sola: el modelo tendrá que combinarlas.
Paso 3: Preparar los datos¶
Dividimos las citas en un conjunto de entrenamiento (75%), del que aprende el modelo, y un conjunto de prueba (25%), que guardamos aparte para evaluarlo con citas que nunca ha visto. stratify=y mantiene la misma proporción de ausencias en los dos conjuntos.
Después escalamos las características: restamos la media y dividimos por la desviación típica, de modo que todas las variables se miden en "desviaciones típicas". Sin esto, una variable en km y otra en días tendrían tamaños muy distintos y el descenso del gradiente convergería despacio. La media y la desviación típica se calculan solo con el conjunto de entrenamiento y se aplican a los dos: el conjunto de prueba no debe influir de ninguna forma en el entrenamiento.
feature_names = ['edad', 'dias_desde_ultima_cita', 'distancia_km', 'hora', 'primera_visita', 'ausencias_previas']
X = df[feature_names].to_numpy(dtype=float)
y = df['no_asiste'].to_numpy()
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, random_state=42, stratify=y)
# Escalado con las estadísticas del conjunto de ENTRENAMIENTO únicamente
mean = X_train.mean(axis=0)
std = X_train.std(axis=0)
X_train_s = (X_train - mean) / std
X_test_s = (X_test - mean) / std
print(f"Entrenamiento: {len(y_train)} citas ({y_train.mean():.1%} de ausencias)")
print(f"Prueba: {len(y_test)} citas ({y_test.mean():.1%} de ausencias)")
Entrenamiento: 900 citas (26.2% de ausencias) Prueba: 300 citas (26.3% de ausencias)
Paso 4: La función sigmoide¶
La sigmoide $\sigma(z) = \frac{1}{1 + e^{-z}}$ convierte la puntuación lineal $z$ en una probabilidad: las puntuaciones muy positivas dan probabilidades cercanas a 1, las muy negativas dan probabilidades cercanas a 0 y $z = 0$ da exactamente 0,5, el umbral por defecto.
def sigmoid(z):
"""Convierte cualquier número real (o array) en un valor entre 0 y 1."""
return 1 / (1 + np.exp(-z))
for z in [-4, -2, 0, 2, 4]:
print(f"sigmoide({z:>2}) = {sigmoid(z):.3f}")
z = np.linspace(-8, 8, 200)
plt.figure(figsize=(8, 4))
plt.plot(z, sigmoid(z), color='tab:orange', linewidth=3)
plt.axhline(0.5, color='gray', linestyle='--', label='Umbral 0,5')
plt.xlabel('Puntuación lineal z = w·x + b')
plt.ylabel('Probabilidad de ausencia')
plt.title('La función sigmoide')
plt.legend()
plt.grid(True)
plt.show()
sigmoide(-4) = 0.018 sigmoide(-2) = 0.119 sigmoide( 0) = 0.500 sigmoide( 2) = 0.881 sigmoide( 4) = 0.982
Paso 5: La función de coste: entropía cruzada (log loss)¶
Para medir cuánto se equivoca el modelo no podemos usar el error cuadrático medio del capítulo 2: con una sigmoide produce una superficie de coste con muchos valles. La regresión logística usa la entropía cruzada, o log loss:
$$J(w, b) = -\frac{1}{m} \sum_{i=1}^{m} \Big[ y_i \log(\hat{p}_i) + (1 - y_i) \log(1 - \hat{p}_i) \Big]$$
Para un paciente que no acudió ($y = 1$) solo cuenta el primer término, $-\log(\hat{p})$: es casi 0 si el modelo le dio una probabilidad cercana a 1 y crece sin límite cuando esa probabilidad se acerca a 0. Los errores cometidos con mucha seguridad se castigan con dureza.
def compute_cost(X, y, w, b):
"""Entropía cruzada (log loss) del modelo con pesos w y sesgo b."""
p = sigmoid(X @ w + b)
p = np.clip(p, 1e-15, 1 - 1e-15) # evita log(0)
return -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))
# El coste de un único paciente que NO acudió (y = 1):
for p in [0.9, 0.5, 0.1]:
print(f"El modelo predijo p = {p}: coste = {-np.log(p):.3f}")
El modelo predijo p = 0.9: coste = 0.105 El modelo predijo p = 0.5: coste = 0.693 El modelo predijo p = 0.1: coste = 2.303
Paso 6: Entrenamiento con descenso del gradiente¶
Los gradientes de la entropía cruzada tienen exactamente la misma forma que en la regresión lineal; la única diferencia es que la predicción pasa ahora por la sigmoide:
$$\frac{\partial J}{\partial w} = \frac{1}{m} X^\top (\hat{p} - y) \qquad \frac{\partial J}{\partial b} = \frac{1}{m} \sum_{i=1}^{m} (\hat{p}_i - y_i)$$
y la regla de actualización es la que ya conocemos: $w := w - \alpha \frac{\partial J}{\partial w}$ y $b := b - \alpha \frac{\partial J}{\partial b}$, donde $\alpha$ es la tasa de aprendizaje. Como las características están escaladas, una tasa de aprendizaje de 0,1 converge rápido.
def train_logistic(X, y, learning_rate=0.1, num_iterations=3000):
"""Descenso del gradiente por lotes sobre la entropía cruzada. Devuelve w, b y el coste de cada iteración."""
m, n = X.shape
w = np.zeros(n)
b = 0.0
costs = []
for i in range(num_iterations):
p = sigmoid(X @ w + b) # a) predecir con los parámetros actuales
dw = (1 / m) * X.T @ (p - y) # b) gradientes
db = (1 / m) * np.sum(p - y)
w = w - learning_rate * dw # c) un paso cuesta abajo
b = b - learning_rate * db
costs.append(compute_cost(X, y, w, b))
if i % 500 == 0:
print(f"Iteración {i:4d}: coste = {costs[-1]:.4f}")
return w, b, costs
w, b, costs = train_logistic(X_train_s, y_train)
print(f"\nCoste final: {costs[-1]:.4f}")
print("Pesos aprendidos (w):", w)
print(f"Sesgo aprendido (b): {b:.3f}")
plt.figure(figsize=(8, 4))
plt.plot(costs)
plt.xlabel('Iteración')
plt.ylabel('Coste (log loss)')
plt.title('Curva de aprendizaje')
plt.grid(True)
plt.show()
Iteración 0: coste = 0.6807 Iteración 500: coste = 0.3523 Iteración 1000: coste = 0.3511 Iteración 1500: coste = 0.3511 Iteración 2000: coste = 0.3511 Iteración 2500: coste = 0.3511 Coste final: 0.3511 Pesos aprendidos (w): [-1.269 1.169 1.02 0.365 0.286 1.026] Sesgo aprendido (b): -1.808
El coste baja rápido y después se aplana: el modelo ha convergido.
Paso 7: Evaluación del modelo¶
Con un umbral convertimos cada probabilidad en una predicción, y al comparar las predicciones con lo que realmente ocurrió hay cuatro resultados posibles, la matriz de confusión:
| Predicción: no acudirá | Predicción: acudirá | |
|---|---|---|
| No acudió | Verdadero positivo (VP) | Falso negativo (FN) |
| Acudió | Falso positivo (FP) | Verdadero negativo (VN) |
A partir de ella calculamos las métricas que se usan en medicina:
- Exactitud (accuracy): $\frac{VP + VN}{\text{total}}$, la proporción de predicciones correctas.
- Sensibilidad (recall): $\frac{VP}{VP + FN}$, la proporción de ausencias que el modelo detecta.
- Especificidad: $\frac{VN}{VN + FP}$, la proporción de asistencias que el modelo reconoce.
- Valor predictivo positivo (VPP): $\frac{VP}{VP + FP}$, cuando el modelo dice "no acudirá", con qué frecuencia acierta.
- Valor predictivo negativo (VPN): $\frac{VN}{VN + FN}$, cuando el modelo dice "acudirá", con qué frecuencia acierta.
Una referencia útil: un "modelo" que predice siempre que todos los pacientes acudirán ya tiene una exactitud igual a la tasa de asistencia, sin detectar ni una sola ausencia. Nuestro modelo tiene que superarlo.
def predict_proba(X, w, b):
return sigmoid(X @ w + b)
def evaluate(y_true, p, threshold=0.5):
"""Matriz de confusión y métricas clínicas para un umbral dado."""
y_pred = (p >= threshold).astype(int)
tp = np.sum((y_pred == 1) & (y_true == 1))
fn = np.sum((y_pred == 0) & (y_true == 1))
fp = np.sum((y_pred == 1) & (y_true == 0))
tn = np.sum((y_pred == 0) & (y_true == 0))
return {
'umbral': threshold, 'VP': tp, 'FN': fn, 'FP': fp, 'VN': tn,
'exactitud': (tp + tn) / len(y_true),
'sensibilidad': tp / (tp + fn),
'especificidad': tn / (tn + fp),
'VPP': tp / (tp + fp),
'VPN': tn / (tn + fn),
}
p_test = predict_proba(X_test_s, w, b)
r = evaluate(y_test, p_test, threshold=0.5)
print("Matriz de confusión (conjunto de prueba, umbral 0,5):")
print(f" Predicción: no acude Predicción: acude")
print(f" No acudió VP = {r['VP']:<17} FN = {r['FN']}")
print(f" Acudió FP = {r['FP']:<17} VN = {r['VN']}")
print()
for k in ['exactitud', 'sensibilidad', 'especificidad', 'VPP', 'VPN']:
print(f"{k:>13}: {r[k]:.3f}")
print(f"\nReferencia (predecir siempre 'acude'): exactitud = {1 - y_test.mean():.3f}")
Matriz de confusión (conjunto de prueba, umbral 0,5):
Predicción: no acude Predicción: acude
No acudió VP = 49 FN = 30
Acudió FP = 22 VN = 199
exactitud: 0.827
sensibilidad: 0.620
especificidad: 0.900
VPP: 0.690
VPN: 0.869
Referencia (predecir siempre 'acude'): exactitud = 0.737
El modelo supera la referencia y es muy fiable cuando dice que un paciente sí acudirá (VPN alto). Pero con el umbral por defecto solo detecta una parte de las ausencias (la sensibilidad): muchas faltas pasan desapercibidas.
Paso 8: Elegir el umbral¶
El umbral de 0,5 es solo una convención. Bajarlo hace que el modelo marque a más pacientes como posibles ausencias: la sensibilidad sube (se escapan menos faltas), pero también aumentan las falsas alarmas, así que bajan la especificidad y el VPP. Qué equilibrio es mejor no es una cuestión matemática, sino clínica y organizativa: aquí una falsa alarma cuesta una llamada de recordatorio, mientras que una ausencia no detectada cuesta un hueco vacío en una lista de espera.
table = pd.DataFrame([evaluate(y_test, p_test, t) for t in [0.5, 0.4, 0.3, 0.2]])
table.round(3)
| umbral | VP | FN | FP | VN | exactitud | sensibilidad | especificidad | VPP | VPN | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.5 | 49 | 30 | 22 | 199 | 0.827 | 0.620 | 0.900 | 0.690 | 0.869 |
| 1 | 0.4 | 54 | 25 | 32 | 189 | 0.810 | 0.684 | 0.855 | 0.628 | 0.883 |
| 2 | 0.3 | 61 | 18 | 46 | 175 | 0.787 | 0.772 | 0.792 | 0.570 | 0.907 |
| 3 | 0.2 | 68 | 11 | 64 | 157 | 0.750 | 0.861 | 0.710 | 0.515 | 0.935 |
Con un umbral de 0,3 el modelo detecta muchas más ausencias, a cambio de más llamadas de recordatorio a pacientes que habrían venido de todos modos.
Paso 9: Interpretar los coeficientes¶
Como las características están escaladas, los pesos son comparables entre sí: cada uno indica cuánto cambian las log-odds (el logaritmo de la razón de probabilidades) de una ausencia cuando esa variable aumenta una desviación típica, con las demás fijas. Su exponencial, $e^{w}$, es la razón de probabilidades (odds ratio): un valor de 2 duplica la probabilidad relativa de ausencia y un valor de 0,5 la reduce a la mitad.
coefficients = pd.DataFrame({'peso': w, 'razón de prob.': np.exp(w)}, index=feature_names)
coefficients = coefficients.reindex(coefficients['peso'].abs().sort_values(ascending=False).index)
print(coefficients.round(3))
plt.figure(figsize=(8, 4))
plt.barh(coefficients.index[::-1], coefficients['peso'][::-1],
color=['tab:red' if v > 0 else 'tab:blue' for v in coefficients['peso'][::-1]])
plt.axvline(0, color='black', linewidth=1)
plt.xlabel('Peso (cambio en log-odds por desviación típica)')
plt.title('Qué hace más probable (rojo) o menos probable (azul) una ausencia')
plt.tight_layout()
plt.show()
peso razón de prob. edad -1.269 0.281 dias_desde_ultima_cita 1.169 3.217 ausencias_previas 1.026 2.791 distancia_km 1.020 2.774 hora 0.365 1.440 primera_visita 0.286 1.332
El modelo ha encontrado el patrón que el equipo sospechaba, y en cabeza aparecen las mismas tres variables que destaca el libro. La edad es la señal más fuerte: cada desviación típica de edad (unos 20 años) multiplica por unas 0,28 la probabilidad relativa de ausencia, así que los pacientes mayores faltan mucho menos. Los días desde la última cita (×3,2), las ausencias previas (×2,8) y la distancia al hospital (×2,8) aumentan el riesgo. La hora de la cita (×1,4) y el tipo de visita (×1,3) importan menos. Son asociaciones aprendidas de los datos, no causas: le dicen al equipo a quién llamar, no por qué un paciente falta a una cita.
Paso 10: Predecir para pacientes nuevos¶
Por último, usamos el modelo entrenado con dos citas nuevas. Los datos nuevos se escalan con la misma media y desviación típica que el conjunto de entrenamiento.
# edad, días desde la última cita, distancia (km), hora, primera visita, ausencias previas
new_patients = np.array([
[35, 100, 18.0, 18, 1, 0], # joven, mucho tiempo sin venir, vive lejos, a última hora, primera visita
[33, 80, 14.0, 17, 1, 0], # parecido, pero algo más cerca del perfil habitual
[62, 30, 6.0, 9, 0, 0], # mayor, revisión reciente, vive cerca, por la mañana
])
p_new = predict_proba((new_patients - mean) / std, w, b)
for patient, p in zip(new_patients, p_new):
flag_05 = "marcado" if p >= 0.5 else "no marcado"
flag_03 = "marcado" if p >= 0.3 else "no marcado"
print(f"Paciente de {patient[0]:.0f} años: P(ausencia) = {p:5.1%} | umbral 0,5: {flag_05:<10} | umbral 0,3: {flag_03}")
Paciente de 35 años: P(ausencia) = 57.5% | umbral 0,5: marcado | umbral 0,3: marcado Paciente de 33 años: P(ausencia) = 36.0% | umbral 0,5: no marcado | umbral 0,3: marcado Paciente de 62 años: P(ausencia) = 0.3% | umbral 0,5: no marcado | umbral 0,3: no marcado
El segundo paciente muestra por qué importa el umbral: con 0,5 el modelo no lo marcaría; con 0,3 el equipo llamaría para confirmar la cita.
Conclusión¶
Hemos construido una regresión logística desde cero: la sigmoide convierte una puntuación lineal en una probabilidad, la entropía cruzada mide el error y el descenso del gradiente ajusta los pesos exactamente igual que en la regresión lineal.
También hemos visto que un clasificador no se evalúa con un solo número. La exactitud oculta que, con el umbral por defecto, muchas ausencias no se detectan; la sensibilidad, la especificidad, el VPP y el VPN cuentan la historia completa, y el umbral permite al equipo elegir el equilibrio que conviene al hospital.
En el siguiente cuaderno, Regresión logística con scikit-learn, entrenaremos el mismo modelo con una biblioteca y añadiremos dos herramientas de este capítulo: la regularización, para evitar el sobreajuste, y la validación cruzada de k particiones, para elegir su intensidad.