pypro.mx

Estructura y componentes de los modelos lineales mixtos

El estudio de sueño tiene 180 filas, pero no son 180 observaciones independientes: son 18 sujetos medidos 10 días cada uno. Dos mediciones del mismo sujeto se parecen más entre sí que dos mediciones de sujetos distintos. Una regresión ordinaria ignora esa estructura y reporta una precisión que no tiene. Un modelo lineal mixto la modela de forma explícita.

Efectos comunes y efectos de grupo

Un modelo mixto combina dos tipos de coeficientes:

Efectos comunes (también conocidos como efectos fijos): parámetros compartidos por toda la población, como el efecto promedio de cada día sin dormir.
Efectos de grupo (específicos de grupo o aleatorios): desviaciones de cada sujeto respecto de ese promedio. No se estiman como parámetros sueltos: se asume que vienen de una distribución común cuya desviación estándar también se estima.

Con interceptos y pendientes aleatorias, para la medición $i$ del sujeto $j$:

$$y_{ij} = \beta_0 + u_{0j} + (\beta_1 + u_{1j})\, x_{ij} + \varepsilon_{ij}$$ $$u_{0j} \sim \mathcal{N}(0, \tau_0), \quad u_{1j} \sim \mathcal{N}(0, \tau_1), \quad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma)$$

Aquí $\beta_0$ y $\beta_1$ son efectos comunes; $u_{0j}$ es el intercepto aleatorio (qué tan lento es el sujeto $j$ en el día 0) y $u_{1j}$ la pendiente aleatoria (cuánto más o menos le afecta cada día). $\tau_0$ y $\tau_1$ miden cuánto varían los sujetos y $\sigma$ es el ruido dentro de cada sujeto. En la versión bayesiana, $\tau_0$ y $\tau_1$ reciben su propio prior: son hiperparámetros.

Tres formas de agrupar

Agrupamiento completo (complete pooling): una sola recta para todos. Ignora las diferencias entre sujetos y mete esa variación en $\sigma$.
Sin agrupamiento (no pooling): una recta independiente por sujeto. Cada estimación usa solo 10 puntos y los sujetos ruidosos dan valores extremos.
Agrupamiento parcial (partial pooling): el modelo mixto. Cada sujeto tiene su recta, pero las estimaciones se jalan hacia el promedio poblacional en proporción a su incertidumbre. Ese jalón se llama shrinkage (encogimiento).

El shrinkage depende de los datos: si $\tau$ es grande frente al ruido de cada grupo, casi no hay encogimiento; si un grupo tiene pocos datos o mucho ruido, se acerca más a la media.

Sintaxis de fórmulas

Bambi usa la sintaxis de lme4 en R, implementada en la librería formulae:

y ~ x + (1|g): intercepto aleatorio por grupo g.
y ~ x + (x|g): intercepto y pendiente de x aleatorios por grupo.
y ~ x + (0 + x|g): solo pendiente aleatoria, sin intercepto de grupo.
y ~ x + (1|g1) + (1|g2): dos agrupamientos cruzados, por ejemplo sujeto y evaluador.
y ~ x + (1|g1/g2): anidados; se expande a (1|g1) + (1|g1:g2), por ejemplo aulas dentro de escuelas.

Aplicación: shrinkage en el estudio de sueño

El siguiente código reconstruye los datos, ajusta el modelo de pendientes aleatorias y compara la pendiente de cada sujeto sin agrupamiento (mínimos cuadrados por sujeto) contra el agrupamiento parcial:

import numpy as np, pandas as pd, bambi as bmb
rng = np.random.default_rng(42)                     # cuerpo de make_sueno()
u0, u1, dias = rng.normal(0, 25, 18), rng.normal(0, 6, 18), np.arange(10)
y = [250 + u0[s] + (10 + u1[s]) * dias + rng.normal(0, 25, 10) for s in range(18)]
df = pd.DataFrame({"sujeto": pd.Categorical(np.repeat([f"S{s+1:02d}" for s in range(18)], 10)),
                   "dias": np.tile(dias, 18), "tiempo_reaccion": np.concatenate(y).round(1)})
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")
post = idata.posterior
parcial = (post["dias"] + post["dias|sujeto"]).mean(("chain", "draw")).to_numpy()
sin_agrup = np.array([np.polyfit(g["dias"], g["tiempo_reaccion"], 1)[0]
                      for _, g in df.groupby("sujeto", observed=True)])
print("S13 y S14 sin agrupar:", sin_agrup[[12, 13]].round(2))
print("S13 y S14 parcial:   ", parcial[[12, 13]].round(2))
print("desv. est. entre sujetos:", sin_agrup.std().round(2), parcial.std().round(2))

El sujeto S13, el más extremo, pasa de 22.58 a 20.63 ms por día; S14 sube de 3.93 a 6.34. La dispersión de las pendientes baja de 5.40 a 4.26. En el posterior, el efecto común dias vale 10.66 (HDI 94%: 8.22 a 13.32), $\tau_1$ (dias|sujeto_sigma) vale 5.23, $\tau_0$ (1|sujeto_sigma) 25.6 y $\sigma$ 22.35. Con agrupamiento completo, dias sale 10.73 con un HDI más angosto (8.80 a 12.50) y $\sigma$ sube a 38.56: el modelo simple esconde la variación entre sujetos en el ruido y sobrestima su certeza sobre el efecto promedio.

Si imprimes modelo verás que Bambi asigna a cada efecto de grupo un prior Normal centrado en 0 cuya desviación tiene a su vez un prior HalfNormal, por ejemplo dias|sujeto ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 42.7091)). Esa escala sale automáticamente de los datos.

Trampas comunes

Usar como grupo una variable con 2 o 3 niveles. Con tan pocos grupos $\tau$ queda mal identificado; suele ser mejor un efecto común.
Agrupar por una variable numérica. El factor de agrupamiento debe ser categórico; convierte los códigos con pd.Categorical.
Olvidar la pendiente aleatoria. Si el efecto varía por grupo y solo modelas el intercepto, el intervalo del efecto común sale artificialmente estrecho.
Leer dias|sujeto como la pendiente del sujeto. Es la desviación; la pendiente del sujeto es dias + dias|sujeto.

Cierre

Antes de escribir la fórmula, identifica la unidad que se repite (sujeto, escuela, tienda) y pregúntate qué coeficientes pueden variar entre unidades. Eso define tus efectos de grupo; el resto son efectos comunes.

Inferencia Bayesiana Aplicada: Modelos Lineales Mixtos con Bambi

Estructura y Componentes de los Modelos