Inferencia Bayesiana Aplicada: Modelos Lineales Mixtos con Bambi

Evaluación de modelos ajustados

Que model.fit() termine sin errores no significa que el posterior sea confiable ni que el modelo describa bien los datos. La evaluación tiene tres capas, en este orden: convergencia (¿las cadenas exploraron bien el posterior?), ajuste (¿el modelo genera datos parecidos a los observados?) y comparación (¿qué modelo predice mejor datos nuevos?). Usaremos el modelo de pendiente aleatoria del estudio de sueño (df = make_sueno(), 18 sujetos × 10 días).

1. Diagnósticos de convergencia

Con MCMC se corren varias cadenas independientes. Si todas describen la misma distribución, el posterior es usable. Las métricas clave aparecen en az.summary:

• R-hat ($\hat R$): compara la varianza entre cadenas con la varianza dentro de cada cadena. Debe estar por debajo de 1.01. Valores como 1.02 o 1.05 indican que las cadenas no coinciden.
• ESS bulk y ESS tail: tamaño de muestra efectivo en el centro y en las colas de la distribución. Las muestras MCMC están autocorrelacionadas, así que 2000 draws valen menos de 2000 muestras independientes. Una regla práctica es pedir al menos 400 (100 por cadena) para confiar en medias e intervalos.
• Divergencias: transiciones donde NUTS no pudo seguir la geometría del posterior. Aunque sean pocas, sesgan las regiones difíciles (típicamente las desviaciones de grupo cercanas a cero).

import bambi as bmb
import arviz as az

kw = dict(draws=1000, tune=1000, chains=2, cores=1,
          random_seed=42, inference_method="nutpie")
m_slope = bmb.Model("tiempo_reaccion ~ dias + (dias|sujeto)", df)
idata_slope = m_slope.fit(**kw)
print(az.summary(idata_slope, var_names=["dias", "sigma",
                 "1|sujeto_sigma", "dias|sujeto_sigma"]))
print("divergencias:", int(idata_slope.sample_stats["diverging"].sum()))
az.plot_trace(idata_slope, var_names=["dias", "dias|sujeto_sigma"]);

Resultado real: dias = 10.66 (HDI 94% de 8.22 a 13.32), sigma = 22.35, 1|sujeto_sigma = 25.59 y dias|sujeto_sigma = 5.23. Todos los R-hat valen 1.00 (el máximo sobre todos los parámetros, incluidos los 36 efectos por sujeto, es 1.005), el ESS bulk mínimo de todo el posterior es 565 (en Intercept) y hay 0 divergencias. En plot_trace busca densidades de ambas cadenas superpuestas y trazas sin tendencias ni tramos atascados.

Si los diagnósticos fallan, las soluciones habituales son más tune y draws, más cadenas, priors más informativos para las desviaciones de grupo o, si hay divergencias, revisar la parametrización. Bambi ya usa por defecto la parametrización no centrada de los efectos de grupo, que evita muchos problemas.

2. Posterior predictive checks

Un posterior bien muestreado puede venir de un modelo equivocado. El chequeo predictivo posterior simula datos nuevos $\tilde y$ desde $p(\tilde y \mid y)$ y los compara con los observados. En Bambi se generan con predict(kind="response"), que agrega el grupo posterior_predictive al InferenceData:

m_slope.predict(idata_slope, kind="response")
az.plot_ppc(idata_slope, num_pp_samples=100);

La gráfica muestra la densidad observada sobre 100 densidades simuladas. Si la observada cae fuera de la nube (colas más pesadas, asimetría, valores imposibles), el modelo no captura algo. Ojo: kind="response_params" (el default) devuelve la media $\mu$ en el grupo posterior, no datos simulados; para un PPC necesitas "response".

3. Comparación con LOO

La validación cruzada leave-one-out estima qué tan bien predice el modelo cada observación sin haberla visto. ArviZ la aproxima con muestreo de importancia suavizado con Pareto (PSIS-LOO) sin reajustar el modelo. La cantidad central es la densidad predictiva log esperada:

$$\text{elpd}_{\text{loo}} = \sum_{i=1}^{n} \log p(y_i \mid y_{-i})$$

Mayor es mejor. Para calcularla hace falta la log-verosimilitud por observación. Gotcha: con inference_method="nutpie" el argumento idata_kwargs={"log_likelihood": True} se ignora; llama a model.compute_log_likelihood(idata) después de ajustar.

formulas = {"pool": "tiempo_reaccion ~ dias",
            "int": "tiempo_reaccion ~ dias + (1|sujeto)",
            "slope": "tiempo_reaccion ~ dias + (dias|sujeto)"}
modelos = {}
for nombre, f in formulas.items():
    m = bmb.Model(f, df)
    idata = m.fit(**kw)
    m.compute_log_likelihood(idata)  # necesario con nutpie
    modelos[nombre] = idata
print(az.compare(modelos))

Resultado real:

• slope: elpd_loo = −829.95, p_loo = 27.4 (rank 0).
• int: elpd_loo = −856.49, p_loo = 18.3, elpd_diff = 26.5, dse = 6.2.
• pool: elpd_loo = −914.06, p_loo = 3.3, elpd_diff = 84.1, dse = 10.1.

Cómo leerlo: elpd_diff es la distancia al mejor modelo y dse su error estándar. Una diferencia de 26.5 con dse de 6.2 son más de 4 errores estándar: la pendiente aleatoria predice claramente mejor. p_loo es el número efectivo de parámetros; aquí 27.4 es mucho menos que los 41 parámetros nominales del modelo, porque el agrupamiento parcial encoge los efectos por sujeto. Si p_loo supera al número real de parámetros o al de observaciones, el modelo está mal especificado.

Por último, revisa Pareto $\hat k$ de cada observación. En ArviZ 0.23 el umbral es $\min(1 - 1/\log_{10} S,\ 0.7)$, con $S$ el número de draws; con 2000 draws equivale a 0.70. Por encima, la aproximación PSIS no es confiable para esa observación y az.compare marca warning=True. En el modelo slope el máximo es 0.63 y ninguna observación supera 0.7.

Trampas comunes

• Mirar solo R-hat. Un R-hat de 1.00 con ESS de 50 da intervalos inestables. Revisa ambos.
• Imprimir az.summary(idata) sin var_names tras compute_log_likelihood. Ese método agrega mu al posterior, una variable por observación, y la tabla se vuelve enorme.
• Elegir por diferencias menores a 2 dse. Ahí los modelos son prácticamente equivalentes; prefiere el más simple o promedia.
• Comparar modelos con distinta respuesta (y contra log(y)): sus elpd no son comparables.

Cierre

Antes de interpretar un solo coeficiente, deja registrados cuatro números: R-hat máximo, ESS mínimo, número de divergencias y el máximo de Pareto $\hat k$. Después agrega la gráfica PPC y la tabla de az.compare. Con eso, la elección del modelo queda justificada con evidencia.