Inferencia Bayesiana Aplicada: Modelos Lineales Mixtos con Bambi

Resumen del curso y conclusiones

En este curso usaste dos conjuntos de datos con estructura de grupos: 18 sujetos medidos durante 10 días de privación de sueño y 1028 alumnos repartidos en 30 escuelas. En ambos casos la pregunta fue la misma: cómo estimar un efecto promedio sin ignorar que las unidades de un mismo grupo se parecen entre sí. Esta lectura reúne el flujo completo, los errores que más se repiten y hacia dónde seguir.

El flujo bayesiano completo

Todo el curso se apoya en el teorema de Bayes, $p(\theta \mid y) \propto p(y \mid \theta)\,p(\theta)$, aplicado a un modelo jerárquico:

$$y_{ij} = (\beta_0 + b_{0j}) + (\beta_1 + b_{1j})\,x_{ij} + \varepsilon_{ij}, \quad b_{j} \sim \mathcal{N}(0, \Sigma_b), \quad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2)$$

Los $\beta$ son los efectos comunes y los $b_j$ los efectos de grupo. Su desviación estándar se estima a partir de los datos, y de ahí sale el agrupamiento parcial: los grupos con pocos datos se acercan a la media y los grupos con muchos datos conservan su propia estimación.

En el estudio de sueño, pasar del modelo pooled al de pendiente aleatoria redujo sigma de 38.56 a 22.35 y mejoró elpd_loo de −914.06 a −829.95. El efecto de dias siguió cerca de 10.7 ms por día, pero su HDI se ensanchó de [8.80, 12.50] a [8.22, 13.32]. Ese intervalo más ancho es más honesto: reconoce que los sujetos responden distinto a la falta de sueño.

Checklist de trabajo con Bambi

1. Explora por grupo: tamaño, media y tendencia de cada grupo antes de modelar.
2. Especifica la fórmula empezando simple: y ~ x, luego y ~ x + (1|g) y luego y ~ x + (x|g).
3. Revisa los priors con print(modelo) y modelo.prior_predictive(). Los priors por defecto de Bambi son débilmente informativos y se escalan a los datos.
4. Ajusta con modelo.fit(). El método por defecto es el NUTS de PyMC; inference_method="nutpie" suele ser más rápido.
5. Diagnostica: R-hat < 1.01, ESS suficiente en bulk y tail, cero divergencias.
6. Valida el ajuste con predict(kind="response") y az.plot_ppc.
7. Compara con compute_log_likelihood y az.compare, leyendo elpd_diff junto con dse y revisando Pareto $\hat k$.
8. Interpreta con HDI, probabilidades posteriores de hipótesis y bmb.interpret, separando la incertidumbre de $\mu$ de la predictiva.

Errores comunes

• Interpretar sin diagnosticar. En el proyecto, el modelo de intercepto aleatorio con 2 cadenas dio R-hat 1.02 y ESS 146 en Intercept. El ajuste terminó sin errores, pero no era utilizable sin más muestras.
• Usar la variable de grupo como efecto común (y ~ x + g) cuando lo que interesa es generalizar a grupos nuevos. Eso estima cada grupo por separado y no permite predecir para uno no visto.
• Olvidar la log-verosimilitud con nutpie. idata_kwargs={"log_likelihood": True} se ignora; sin compute_log_likelihood, az.compare falla.
• Dejar que bmb.interpret elija el grupo. Si no indicas el grupo en conditional, Bambi fija uno por defecto y las gráficas no muestran lo que crees.
• Confundir intervalos. El HDI de $\mu$ describe la media; para una observación nueva se usa el predictivo, que es más ancho.

Cuándo usar modelos mixtos

Úsalos cuando los datos tienen grupos (sujetos, escuelas, tiendas, pacientes, regiones), medidas repetidas o grupos de tamaño desigual. Son especialmente útiles cuando algunos grupos tienen pocos datos o cuando quieres predecir para grupos nuevos. Si solo hay 2 o 3 grupos, la desviación entre grupos queda mal identificada y un efecto común con prior informativo suele bastar.

Extensiones

• GLMM: la misma sintaxis sirve para respuestas no gaussianas. family="bernoulli" para resultados binarios y family="poisson" para conteos, con los mismos términos (x|g).
• Modelos distribucionales: con bmb.Formula("y ~ x + (1|g)", "sigma ~ x") también modelas la varianza, útil cuando la dispersión crece con el predictor.
• Splines: bs(x, df=4) dentro de la fórmula agrega relaciones no lineales sin elegir una forma funcional.
• Procesos gaussianos aproximados: hsgp(x, m=10, c=2) agrega una función suave flexible con un costo computacional manejable.

Las cuatro extensiones se construyen con bmb.Model en Bambi 0.17.2 sin código adicional de PyMC.

Siguientes pasos de lectura

• Andrew Gelman y Jennifer Hill, Data Analysis Using Regression and Multilevel/Hierarchical Models (2007): la referencia clásica de modelos multinivel.
• Richard McElreath, Statistical Rethinking, 2.ª edición (2020): una introducción conceptual a modelos bayesianos y jerárquicos.
• Capretto et al., "Bambi: A Simple Interface for Fitting Bayesian Linear Models in Python", Journal of Statistical Software 103(15), 2022.
• La documentación de Bambi, con ejemplos de GLMM, modelos distribucionales, splines y HSGP.
• La documentación de ArviZ para diagnósticos y LOO.

Cierre

La idea central del curso: un buen análisis bayesiano no termina en el ajuste. Termina cuando puedes defender el modelo con diagnósticos, chequeos predictivos y una comparación formal, y cuando cada conclusión lleva su intervalo de credibilidad.