Inferencia Bayesiana Aplicada: Modelos Lineales Mixtos con Bambi

Ejemplo práctico de especificación

En Bambi, especificar un modelo es escribir una fórmula, y esa línea produce un modelo bayesiano completo: verosimilitud, enlace y un prior por parámetro. En este notebook especificas los tres modelos del curso con el estudio de sueño, revisas los priors por defecto, los reemplazas por priors propios y los pones a prueba con una simulación previa a los datos. Descarga el notebook adjunto y ejecútalo celda por celda.

Los datos

El notebook genera, sin descargar nada, 18 sujetos (S01 a S18) medidos durante 10 días de privación de sueño (días 0 a 9): 180 filas con sujeto, dias y tiempo_reaccion en milisegundos. El tiempo de reacción medio es 297.5 ms con desviación estándar 49.2. La columna sujeto es de tipo Categorical; eso le indica a Bambi que es un factor de agrupamiento y no un número.

Paso 1: tres fórmulas, tres supuestos

La sintaxis sigue la notación de lme4 en R. El término (1|sujeto) añade un intercepto por sujeto; (dias|sujeto) añade intercepto y pendiente por sujeto:

import bambi as bmb

pooled = bmb.Model("tiempo_reaccion ~ dias", df)
interc = bmb.Model("tiempo_reaccion ~ dias + (1|sujeto)", df)
slope = bmb.Model("tiempo_reaccion ~ dias + (dias|sujeto)", df)
print(slope)

En notación matemática, el modelo de pendiente aleatoria es

$$y_{ij} = (\beta_0 + u_{0j}) + (\beta_1 + u_{1j})\,\text{dias}_{ij} + \varepsilon_{ij}, \quad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma)$$

donde $\beta_0, \beta_1$ son los efectos comunes y $u_{0j}, u_{1j}$ los efectos de grupo del sujeto $j$. El modelo pooled ignora a los sujetos; interc permite que cada uno parta de un nivel distinto, y slope permite además que cada uno se deteriore a su propio ritmo. Otras variantes útiles: (0 + dias|sujeto) da solo pendientes por sujeto, y 0 + fuera del paréntesis elimina el intercepto común.

Paso 2: leer los priors por defecto

print(slope) muestra, entre otras líneas:

Intercept ~ Normal(mu: 297.5078, sigma: 228.0042)
dias ~ Normal(mu: 0.0, sigma: 42.7091)
1|sujeto ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 228.0042))
dias|sujeto ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 42.7091))
sigma ~ HalfStudentT(nu: 4.0, sigma: 49.069)

Bambi construye priors débilmente informativos escalados a los datos: el intercepto se centra en la media de la respuesta (297.5) y la escala del prior de sigma (49.07) es la desviación estándar de la respuesta. Los efectos de grupo son normales con media cero y una desviación estándar que también se estima, con su propio prior HalfNormal. Esa estructura jerárquica es la que produce el agrupamiento parcial.

Paso 3: categóricas, familia y enlace

Una variable de texto o Categorical se codifica con $N-1$ indicadoras; el primer nivel queda como referencia dentro del intercepto. Si creas fase con los niveles adaptacion y privacion, el coeficiente estimado es el de privacion. Para tratar un número como factor se usa C(dias) en la fórmula o el argumento categorical=["dias"].

La familia define la distribución de la respuesta y, si no indicas link, Bambi usa el enlace canónico. Con una respuesta binaria o un conteo, la fórmula no cambia:

df["lento"] = (df["tiempo_reaccion"] > 300).astype(int)
m_bin = bmb.Model("lento ~ dias + (1|sujeto)", df, family="bernoulli")
print(m_bin)  # Link: p = logit
# para conteos: family="poisson" usa Link: mu = log

Paso 4: priors personalizados con bmb.Prior

Supón que estudios previos sugieren un deterioro cercano a 10 ms por día. Se pasa un diccionario cuyo nombre de término coincide con el que imprime el modelo:

priors = {
    "dias": bmb.Prior("Normal", mu=10, sigma=5),
    "sigma": bmb.Prior("HalfNormal", sigma=30),
    "1|sujeto": bmb.Prior("Normal", mu=0,
                          sigma=bmb.Prior("HalfNormal", sigma=20)),
}
custom = bmb.Model("tiempo_reaccion ~ dias + (1|sujeto)", df,
                   priors=priors)
print(custom)

El prior de un efecto de grupo es anidado: la desviación estándar entre sujetos lleva su propio bmb.Prior. Observa un efecto lateral: al fijar el prior de dias, el prior por defecto del intercepto cambió su sigma de 228.0 a 122.7, porque Bambi lo calcula a partir de los demás términos. Revisa siempre print después de modificar priors.

Paso 5: build() y predictiva previa

model.build() crea el modelo de PyMC sin muestrear. Es obligatorio antes de prior_predictive; sin él, Bambi lanza ValueError: Model is not built yet!.

import numpy as np

custom.build()
pp = custom.prior_predictive(draws=500, random_seed=42)
y_prior = pp.prior_predictive["tiempo_reaccion"].values
print(np.percentile(y_prior, [3, 50, 97]).round(1))
print("fracción negativa:", (y_prior < 0).mean().round(3))

Con los priors personalizados, el 94% central de los tiempos simulados va de 60.0 a 548.1 ms y solo el 1% es negativo. Con los priors por defecto el mismo intervalo va de −353.8 a 976.5 ms y el 17% de los valores es negativo: los priors por defecto son deliberadamente amplios.

Al ajustar con nutpie, los tres modelos reproducen los resultados del curso: la pendiente de dias es 10.73 (pooled), 10.76 (intercepto aleatorio) y 10.66 (pendiente aleatoria), con $\sigma$ bajando de 38.56 a 26.74 y a 22.35. El modelo personalizado da 10.73 para dias, casi igual al de priors por defecto: con 180 observaciones, el prior pesa poco.

Trampas comunes

• Dejar sujeto como número entero: Bambi lo trataría como predictor continuo.
• Escribir mal la clave del diccionario de priors: debe coincidir con el término, por ejemplo "1|sujeto".
• Llamar a prior_predictive sin build() previo.
• Confundir (dias|sujeto) con (0 + dias|sujeto): el primero incluye interceptos por sujeto.

Cierre

Antes de ajustar cualquier modelo, sigue tres pasos: print(model), build() y una predictiva previa. Si los datos simulados son físicamente imposibles a gran escala, ajusta los priors antes de mirar el posterior.