Regresión Logística Bayesiana con PyMC

Ejemplo práctico de regresión logística

En este notebook recorres el flujo completo de una regresión logística bayesiana sobre las 800 solicitudes de crédito del curso: ajustar el modelo, revisar que el muestreo sea confiable, traducir los coeficientes a un lenguaje de negocio y predecir el riesgo de un cliente nuevo con su incertidumbre. Descarga el notebook adjunto y ejecútalo celda por celda.

1. Ajustar el modelo

Partimos del modelo de la lección anterior, con los predictores estandarizados guardados en un pm.Data para poder reemplazarlos después:

cols = ["ingreso", "deuda", "antiguedad"]
media, desv = df[cols].mean(), df[cols].std()
Xs = ((df[cols] - media) / desv).values

with pm.Model(coords={"pred": cols}) as modelo:
    X = pm.Data("X", Xs)
    alpha = pm.Normal("alpha", 0, 1.5)
    beta = pm.Normal("beta", 0, 1, dims="pred")
    pm.Bernoulli("impago", logit_p=alpha + X @ beta,
                 observed=df["impago"].values, shape=X.shape[0])
    idata = pm.sample(1000, tune=1000, chains=4, cores=1, random_seed=42)

El argumento shape=X.shape[0] hace que el tamaño de la verosimilitud siga al de X: cuando cambies los datos por un solo cliente, el modelo generará una sola predicción.

2. Diagnóstico del muestreo

Antes de interpretar nada, verifica tres cosas en az.summary(idata) y en idata.sample_stats:

• Divergencias: 0. Una divergencia indica que NUTS no pudo explorar bien alguna región de la posterior.
• $\hat{R}$: todas ≤ 1.002. Valores por encima de 1.01 indican que las cadenas no coinciden entre sí.
• ESS: más de 4,300 muestras efectivas por parámetro, de 4,000 sorteos. Con 400 o más ya es suficiente para resúmenes estables.

az.plot_trace completa la revisión: las cuatro cadenas deben verse como «orugas peludas» superpuestas, sin tendencias.

3. Interpretar los coeficientes

En escala estandarizada la posterior da: ingreso −0.651 (HDI 94 % [−0.817, −0.461]), deuda 0.625 ([0.463, 0.786]) y antigüedad −0.711 ([−0.878, −0.537]). Para hablar en unidades reales se divide cada muestra entre la desviación estándar del predictor y se exponencia:

beta_orig = idata.posterior["beta"] / desv.values
razon_momios = np.exp(beta_orig.mean(("chain", "draw")))
print(beta_orig.mean(("chain", "draw")).values.round(4))
print(razon_momios.values.round(3))

Los coeficientes en escala original son −0.0826, 3.857 y −0.1492, muy cerca de los valores con los que se generaron los datos (−0.08, 4 y −0.15). En razones de momios: cada mil pesos más de ingreso multiplica los momios de impago por 0.921 (8 % menos); cada año de antigüedad, por 0.862; y subir la razón deuda/ingreso 0.1 los multiplica por $e^{0.3857} \approx 1.47$. Como la posterior son muestras, también puedes afirmar que $P(\beta_{\text{deuda}} > 0 \mid \text{datos}) = 1.0$ en las 4,000 muestras.

4. Predecir un cliente nuevo

Un solicitante gana 18 mil pesos, tiene una razón deuda/ingreso de 0.45 y lleva 2 años en su empleo. Lo estandarizamos con la media y la desviación de entrenamiento, reemplazamos los datos y generamos predicciones:

nuevo = (np.array([[18.0, 0.45, 2]]) - media.values) / desv.values

with modelo:
    pm.set_data({"X": nuevo})
    pred = pm.sample_posterior_predictive(idata, predictions=True, random_seed=42)

eta = idata.posterior["alpha"] + (idata.posterior["beta"] * nuevo[0]).sum("pred")
p_nuevo = expit(eta.values.ravel())
print(p_nuevo.mean().round(3), az.hdi(p_nuevo, hdi_prob=0.94).round(3))

La probabilidad de impago tiene media 0.717 y HDI del 94 % entre 0.652 y 0.782. El 72.1 % de los resultados simulados con predictions=True son impagos, lo que coincide con la probabilidad media. La diferencia frente a un modelo clásico es el intervalo: sabes no solo que el riesgo es alto, sino cuánta incertidumbre tiene esa estimación.

5. Comprobación predictiva posterior

Por último, pm.sample_posterior_predictive(idata, extend_inferencedata=True) simula conjuntos de datos completos con la posterior. La tasa de impago simulada promedio es 0.298, frente a 0.2975 observada. El notebook lo dibuja con az.plot_ppc.

Trampas comunes

• Interpretar sin diagnosticar. Un resumen con divergencias o $\hat{R}$ alto no es confiable, aunque sus números se vean razonables.
• Olvidar shape=X.shape[0]. Sin él, al cambiar X por un solo cliente la verosimilitud conserva sus 800 posiciones y la predicción falla por formas incompatibles.
• Exponenciar la media en lugar de las muestras cuando necesitas un intervalo: calcula el HDI de np.exp(beta) muestra por muestra.
• Reportar solo la media de la predicción. El valor del enfoque bayesiano está en el intervalo; repórtalo siempre.

Cierre

Ajustar, diagnosticar, interpretar y predecir con incertidumbre: ese ciclo se repite en cualquier regresión logística bayesiana, y es la base del proyecto final.