Inferencia Bayesiana Aplicada: Modelos Lineales Mixtos con Bambi

Incertidumbre y credibilidad

El resultado de un modelo bayesiano no es un número sino una distribución. Este notebook muestra cómo convertir ese posterior en respuestas concretas: intervalos de credibilidad, probabilidades de hipótesis, predicciones para sujetos conocidos y nuevos, y efectos por sujeto. Seguimos con el estudio de sueño y el modelo de pendiente aleatoria, que en la lección anterior fue el mejor por LOO.

1. Datos y modelo

El generador es el mismo de todo el curso, así que el notebook es autocontenido:

import numpy as np, pandas as pd, bambi as bmb, arviz as az

def make_sueno(seed=42):
    rng = np.random.default_rng(seed)
    u0, u1 = rng.normal(0, 25, 18), rng.normal(0, 6, 18)
    dias, filas = np.arange(10), []
    for s in range(18):
        y = 250 + u0[s] + (10 + u1[s]) * dias + rng.normal(0, 25, 10)
        filas += [{"sujeto": f"S{s+1:02d}", "dias": int(d), "tiempo_reaccion": round(float(v), 1)}
                  for d, v in zip(dias, y)]
    df = pd.DataFrame(filas)
    df["sujeto"] = pd.Categorical(df["sujeto"])
    return df

df = make_sueno()
modelo = bmb.Model("tiempo_reaccion ~ dias + (dias|sujeto)", df)
idata = modelo.fit(draws=1000, tune=1000, chains=2, cores=1,
                   random_seed=42, inference_method="nutpie")

2. HDI contra ETI y probabilidad de una hipótesis

Hay dos formas comunes de resumir un posterior con un intervalo. El HDI (intervalo de máxima densidad) es el intervalo más corto que contiene la probabilidad pedida; es lo que reporta az.summary (columnas hdi_3% y hdi_97%). El ETI (intervalo de colas iguales) deja la misma probabilidad fuera de cada lado, es decir, usa los percentiles 3 y 97. En distribuciones simétricas casi coinciden; en distribuciones asimétricas, como las desviaciones estándar, se separan.

b = idata.posterior["dias"].values.ravel()
print("HDI 94%:", az.hdi(b, hdi_prob=0.94).round(2))
print("ETI 94%:", np.percentile(b, [3, 97]).round(2))
s = idata.posterior["dias|sujeto_sigma"].values.ravel()
print("sd pendiente  HDI:", az.hdi(s, hdi_prob=0.94).round(2),
      " ETI:", np.percentile(s, [3, 97]).round(2))
for umbral in [8, 10, 12]:
    print(f"P(dias > {umbral}) = {(b > umbral).mean():.3f}")

Para el efecto común dias, el HDI es [8.22, 13.32] y el ETI [8.04, 13.21]: casi iguales. Para dias|sujeto_sigma, el HDI es [3.31, 7.53] y el ETI [3.44, 7.84]: la cola derecha larga desplaza el ETI hacia arriba. La interpretación sí es directa: dados el modelo y los datos, hay 94% de probabilidad de que el efecto promedio esté en ese intervalo.

Una hipótesis como "cada día sin dormir agrega más de 8 ms" se responde contando draws: $P(\beta_{dias} > 8 \mid y) = 0.970$. Para umbrales de 10 y 12 ms resulta 0.694 y 0.166. Es una probabilidad sobre el parámetro, algo que un valor p no ofrece.

3. Efectos de grupo y shrinkage

La pendiente de cada sujeto es $\beta_{dias} + b_{dias,j}$. Si la comparas con una regresión independiente por sujeto (OLS), verás el agrupamiento parcial en acción:

post = idata.posterior
b_mix = (post["dias"] + post["dias|sujeto"]).mean(("chain", "draw")).to_series()
b_ols = df.groupby("sujeto", observed=True).apply(
    lambda g: np.polyfit(g["dias"], g["tiempo_reaccion"], 1)[0])
tabla = pd.DataFrame({"ols": b_ols, "mixto": b_mix.values}).round(2)
print(tabla.loc[["S05", "S11", "S13", "S14", "S18"]])
print("sd OLS:", round(b_ols.std(), 2), " sd mixto:", round(b_mix.std(), 2))

Las pendientes extremas se acercan a la media: S13 pasa de 22.58 a 20.63, S18 de 20.36 a 17.24 y S14 de 3.93 a 6.34. La dispersión entre sujetos baja de 5.56 a 4.38. Con solo 10 puntos por sujeto, una pendiente OLS extrema se debe en parte al ruido; el modelo la corrige usando la información de los otros 17 sujetos.

4. Incertidumbre de parámetros contra incertidumbre predictiva

Hay dos preguntas distintas: "¿cuál es el tiempo medio de S01 en el día 9?" (kind="response_params", la media $\mu$) y "¿cuánto medirá la próxima prueba de S01?" (kind="response", que suma el ruido $\sigma$). Además, un sujeto nuevo (S99) solo se puede predecir con sample_new_groups=True.

nuevo = pd.DataFrame({"sujeto": ["S01", "S99"], "dias": [9, 9]})
opts = dict(data=nuevo, inplace=False, sample_new_groups=True, random_seed=42)
mu = modelo.predict(idata, kind="response_params", **opts).posterior["mu"]
y = modelo.predict(idata, kind="response", **opts).posterior_predictive["tiempo_reaccion"]
mu, y = mu.values.reshape(-1, 2), y.values.reshape(-1, 2)
for j, suj in enumerate(["S01", "S99"]):
    print(suj, "media:", mu[:, j].mean().round(1),
          "HDI mu:", az.hdi(mu[:, j], hdi_prob=0.94).round(1),
          "HDI y:", az.hdi(y[:, j], hdi_prob=0.94).round(1))

Para S01: media 397.7 ms, HDI de $\mu$ [374.5, 419.4] y HDI predictivo [353.1, 449.0]. Para el sujeto nuevo: media 346.8, HDI de $\mu$ [258.9, 412.9] y predictivo [247.6, 431.7]. El HDI de $\mu$ de S99 es unas 3.4 veces más ancho que el de S01 porque no conocemos su intercepto ni su pendiente. En Bambi 0.17, cada draw de un grupo nuevo toma los efectos de un grupo existente elegido al azar (equivale a sample_new_levels="uncertainty" de brms). Si solo quieres la curva poblacional, usa include_group_specific=False: da 345.1 con HDI [321.0, 369.5].

5. Predicciones y contrastes con bmb.interpret

comp = bmb.interpret.comparisons(modelo, idata, contrast={"dias": [0, 9]},
                                 conditional={"sujeto": ["S06", "S13"]})
print(comp[["sujeto", "estimate", "lower_3.0%", "upper_97.0%"]].round(1))
bmb.interpret.plot_predictions(modelo, idata,
    conditional={"dias": np.arange(10), "sujeto": ["S06", "S13", "S18"]});

Nueve días sin dormir suman 48.2 ms a S06 (HDI [15.3, 84.6]) y 185.7 ms a S13 ([149.2, 223.5]). La gráfica dibuja la recta media de cada sujeto con su banda de credibilidad.

Trampas comunes

• Omitir el grupo en bmb.interpret. Con plot_predictions(modelo, idata, ["dias"]), Bambi fija sujeto en un valor por defecto (S01) y solo lo avisa con el mensaje "Default computed for unspecified variable". Indica los sujetos en conditional.
• Reportar el intervalo de $\mu$ como si fuera de una observación. El predictivo siempre es más ancho.
• Leer el HDI como un intervalo de confianza. Es una afirmación de probabilidad sobre el parámetro, condicionada al modelo y a los priors.