Regresión Logística Bayesiana con PyMC

Tu primer modelo logístico: priors y predictiva a priori

Antes de ajustar un modelo a los datos conviene preguntarle qué cree que va a pasar. En este notebook construyes la regresión logística bayesiana desde sus piezas (la función logística, la verosimilitud Bernoulli y los priors) y haces una comprobación predictiva a priori: simulas datos con el modelo sin haber visto ninguna observación y verificas que las predicciones sean razonables. Descarga el notebook adjunto y ejecútalo celda por celda.

De una recta a una probabilidad

La regresión lineal predice cualquier número real; una probabilidad tiene que estar entre 0 y 1. La regresión logística resuelve esto en dos pasos. Primero calcula un predictor lineal $\eta = \alpha + \beta_1 x_1 + \dots + \beta_k x_k$, que vive en la escala de los log-odds. Después lo pasa por la función logística (sigmoide):

$$p = \sigma(\eta) = \frac{1}{1 + e^{-\eta}}, \qquad \operatorname{logit}(p) = \log\frac{p}{1-p} = \eta$$

Cada observación $y_i \in \{0, 1\}$ sigue una Bernoulli con probabilidad $p_i$. Un coeficiente $\beta_j$ dice cuánto cambian los log-odds cuando $x_j$ sube una unidad; $e^{\beta_j}$ es la razón de momios (odds ratio).

Los datos del curso

El notebook genera 800 solicitudes de crédito con semilla fija. Como conocemos los coeficientes verdaderos, podremos comprobar si la posterior los recupera:

import numpy as np
import pandas as pd
from scipy.special import expit

def generar_creditos(n=800, semilla=12):
    rng = np.random.default_rng(semilla)
    ingreso = rng.normal(25, 8, n).clip(5, None)   # miles de pesos al mes
    deuda = rng.beta(2, 5, n)                       # razón deuda/ingreso
    antiguedad = rng.integers(0, 16, n)             # años en el empleo
    z = -1.0 - 0.08 * (ingreso - 25) + 4.0 * (deuda - 0.3) - 0.15 * (antiguedad - 7)
    impago = rng.binomial(1, expit(z))
    return pd.DataFrame(dict(ingreso=ingreso.round(1), deuda=deuda.round(3),
                             antiguedad=antiguedad, impago=impago))

La tasa de impago observada es 0.2975. Antes de modelar estandarizamos los tres predictores (restar la media y dividir entre la desviación estándar), así un prior $\mathcal{N}(0, 1)$ tiene el mismo significado para todos.

Priors con sentido

Un prior $\beta \sim \mathcal{N}(0, 1)$ sobre un predictor estandarizado dice: «es poco probable que subir una desviación estándar multiplique los momios por más de $e^2 \approx 7.4$». Para el intercepto usamos $\alpha \sim \mathcal{N}(0, 1.5)$, que deja la tasa base libre entre casi 0 y casi 1. Es tentador usar priors «no informativos» muy anchos, como $\mathcal{N}(0, 10)$, pero en escala logit eso no es neutral: empuja las probabilidades hacia los extremos.

Comprobación predictiva a priori

pm.sample_prior_predictive sortea parámetros del prior y simula datos con ellos:

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")
    p = pm.Deterministic("p", pm.math.invlogit(alpha + X @ beta))
    pm.Bernoulli("impago", p=p, observed=y)
    prior = pm.sample_prior_predictive(1000, random_seed=42)

pp = prior.prior["p"].values.ravel()
print(np.mean((pp < 0.01) | (pp > 0.99)))   # fracción de probabilidades extremas

Con los priors del curso, solo el 4.5 % de las probabilidades simuladas son extremas (menores que 0.01 o mayores que 0.99) y el 67.3 % cae entre 0.1 y 0.9; la tasa de impago simulada promedio es 0.505. Con $\mathcal{N}(0, 10)$ en todos los parámetros, el 80.2 % de las probabilidades son extremas y solo el 9.6 % queda entre 0.1 y 0.9: el modelo «cree» que casi todos los clientes pagan seguro o no pagan seguro. El notebook dibuja ambos histogramas para que veas la diferencia.

Primera posterior

Con el mismo modelo, pm.sample(1000, tune=1000, chains=4, random_seed=42) tarda unos 4 segundos. La posterior queda cerca de los valores verdaderos en escala estandarizada: ingreso −0.651 (verdad −0.631), deuda 0.625 (verdad 0.648) y antigüedad −0.711 (verdad −0.715). El notebook compara prior y posterior de cada coeficiente en la misma gráfica: el prior es ancho y la posterior, mucho más estrecha, porque 800 observaciones aportan mucha información.

Trampas comunes

• Creer que un prior ancho es neutral. En escala logit, la amplitud se traduce en probabilidades extremas; revísalo con la predictiva a priori.
• Aplicar la sigmoide y usar logit_p a la vez. Usa p=invlogit(eta) o logit_p=eta, nunca las dos transformaciones.
• Estandarizar con estadísticas del conjunto completo y luego predecir sin guardar la media y la desviación: las necesitarás para transformar los clientes nuevos.

Cierre

La predictiva a priori es barata y evita modelos que «saben» cosas absurdas antes de ver los datos. Hazla siempre antes de muestrear la posterior.