8  Estimación de Modelos Bayesianos (Parte 2)

8.1 ANOVA de una vía

library(tidyverse)
library(ggformula)
library(dplyr)
library(R2jags)
library(CalvinBayes)
library(bayesplot)
library(mosaic)
library(rstan)
library(mosaicData)
library(brms)
options(width = 100)

La especie estudiada fue la mosca de la fruta, Drosophila melanogaster (Hanley & Shapiro, 1994), conocida porque las hembras recién inseminadas no se aparean nuevamente durante al menos dos días y los machos no cortejan activamente a las hembras embarazadas. Para ilustrar el uso del modelo, se estudió la duración de la vida de los machos en función de su actividad sexual. En el experimento, se manipuló la actividad sexual suministrando a los machos individuales hembras vírgenes receptivas a una tasa de una u ocho vírgenes por día. La longevidad de estos machos se registró y se comparó con la de dos grupos de control: uno con machos mantenidos con hembras recién inseminadas en igual número que las hembras vírgenes y otro con machos sin hembras. El objetivo era determinar si la actividad sexual de los machos reducía su vida, ya que esto ya estaba establecido para las hembras. Un efecto perjudicial de la actividad sexual en machos sería sorprendente, dado que se presume que los costos fisiológicos de la actividad sexual son menores en machos que en hembras. Cada grupo experimental y de control constaba de 25 moscas macho.

gf_violin(longevity ~ group, data = FruitflyReduced) %>%
  gf_jitter(width = 0.2, height = 0, alpha = 0.5) %>%
  gf_point(stat = "summary", color = "red", size = 3, alpha = 0.5, fun = mean)

8.1.1 Modelo 1

Es bastante fácil pedirle a brm() que ajuste un modelo para nosotros. Solo proporcionemos nuestras variables explicativas y de respuesta y veamos qué sucede.

flies_brm <- brm(longevity ~ group, data = FruitflyReduced)

flies_stan <- stanfit(flies_brm)
flies_stan
mcmc_combo(as.mcmc.list(flies_stan))
mcmc_areasridges(as.mcmc.list(flies_stan), regex_pars = "b_g")

Analicemos la codificación que usa STAN en la variable categórica:

standata(flies_brm) %>% lapply(head)

El modelo resulta ser en este caso:

\[ \mathrm{longevity} = \beta_0 \cdot 1 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \mathrm{noise} \\ = \beta_0 \cdot 1 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + {\sf Norm}(0, \sigma) \]

donde, por ejemplo,

\[ x_1 \;=\; \begin{cases} 1, & \text{si } \text{group} = \text{Pregnant1},\\ 0, & \text{si } \text{group} \neq \text{Pregnant1}. \end{cases} \]

En otras palabras, la distribución de la longevidad es

• \({\sf Norm}(\beta_0, \sigma)\) para el grupo None0
• \({\sf Norm}(\beta_0 + \beta_1, \sigma)\) para el grupo Pregnant1
• \({\sf Norm}(\beta_0 + \beta_2, \sigma)\) para el grupo Pregnant2
• \({\sf Norm}(\beta_0 + \beta_3, \sigma)\) para el grupo Virgin1
• \({\sf Norm}(\beta_0 + \beta_4, \sigma)\) para el grupo Virgin2

Aquí están las distribuciones previas predeterminadas:

prior_summary(flies_brm)

• Previas planas e impropias para los parámetros b_.

• Distribución t con 3 grados de libertad para el intercepto (colas más pesadas que una distribución normal).

  • Nota: Esto es realmente una previa para \(\alpha_0\) (efecto medio de toda la población), no \(\beta_0\), ya que usualmente es más fácil especificar una previa para \(\alpha_0\). Si por alguna razón quisiéramos especificar una previa para \(\beta_0\) en su lugar, hay un pequeño truco: usar la fórmula longevity ~ 0 + intercept + group.
  • La procedencia de los valores 58 y 18 se pueden explicar en lo siguiente:

  df_stats(~ longevity, data = FruitflyReduced, mean, sd)

• “T” para sigma.

8.1.2 Modelo 2: Previas personalizadas

El modelo anterior difiere del modelo básico en el Kruschke.

Podemos acercarnos más al modelo de DBDA usando esto:

flies2_brm <- 
  brm(longevity ~ group, data = FruitflyReduced,
      prior = c(
        set_prior(class = "Intercept", "normal(60, 100)"),  # 100 = 5 * 20
        set_prior(class = "b", "normal(0, 10)"),  # group = "b" es el valor predeterminado; se podría omitir
        set_prior(class = "sigma", "uniform(20.0/1000.0, 20.0 * 1000.0)")
      )
  )
prior_summary(flies2_brm)
stancode(flies2_brm)

Esto todavía no es exactamente igual al modelo utilizado por Kruschke. Resulta que hay múltiples maneras de codificar los \(\beta\)s. Un tercero es utilizado por Kruschke y requiere un poco más de trabajo para ajustarse usando brm():

8.1.3 Modelo 3

Si eliminamos el intercepto en el modelo brm(), obtenemos un modelo con un \(\beta_i\) para cada media de grupo en lugar de un \(\beta_0\) para el primer grupo y un \(\beta_i\) para la diferencia en las medias de los grupos cuando \(i > 0\):

flies3_brm <- 
  brm(
    longevity ~ 0 + group, data = FruitflyReduced,
    prior = c(
      set_prior(class = "b", "normal(60, 10)"),  # group = "b" es el valor predeterminado; se podría omitir
      set_prior(class = "sigma", "uniform(20.0/1000.0, 20.0 * 1000.0)")
    ),
    sample_prior = TRUE
  )
prior_summary(flies3_brm)
stancode(flies3_brm)

Esto es equivalente a

\[\begin{align*} \mathrm{longevity} = \beta_0 x_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \mathrm{ruido} = \beta_0 x_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + {\sf Norm}(0, \sigma) \end{align*}\]

donde, por ejemplo,

\[ x_1 = \begin{cases} 1, & \text{si } \text{group} = \text{Pregnant1},\\ 0, & \text{si } \text{group} \neq \text{Pregnant1}. \end{cases} \]

En otras palabras, la distribución de la longevidad es

• \({\sf Norm}(\beta_0, \sigma)\) para el grupo None0
• \({\sf Norm}(\beta_1, \sigma)\) para el grupo Pregnant1
• \({\sf Norm}(\beta_2, \sigma)\) para el grupo Pregnant2
• \({\sf Norm}(\beta_3, \sigma)\) para el grupo Virgin1
• \({\sf Norm}(\beta_4, \sigma)\) para el grupo Virgin2

8.1.4 Comparando grupos

8.1.4.1 Comparación con el “grupo de referencia”

Usamos el modelo 2:

stanfit(flies2_brm)

En este modelo, un grupo corresponde al intercepto del modelo, y las comparaciones de otros grupos con este grupo implican investigar la distribución posterior de uno de los otros \(\beta\)s.

flies_post <- posterior(flies_stan)
names(flies_post)

plot_post(flies_post$b_groupPregnant1)
hdi(flies_post$b_groupPregnant1)
plot_post(flies_post$b_groupVirgin1)
mcmc_areas(as.mcmc.list(flies_stan), pars = "b_groupVirgin1", prob = 0.95)
hdi(flies_post$b_groupVirgin1)

8.1.4.2 Comparación de otros pares de grupos

¿Qué pasa si queremos comparar los grupos Virgin1 y Virgin8? Podemos usar la identidad \(\beta_0 + \beta_i - (\beta_0 + \beta_j) = \beta_i - \beta_j\) para simplificar el álgebra y hacerlo de esta manera.

flies_post <-
  flies_post %>% mutate(dVirgin = b_groupVirgin8 - b_groupVirgin1)
plot_post(flies_post$dVirgin, xlab = "Virgin8 - Virgin1")

8.1.4.3 Contrastes: Comparación de “grupos de grupos”

¿Qué pasa si queremos comparar los dos grupos de vírgenes con los otros 3 grupos? Esto es un poco más fácil de hacer usando el modelo sin un término de intercepto.

flies3_post <- posterior(flies3_brm)
names(flies3_post)
flies3_post <-
  flies3_post %>% 
  mutate(
    contrast = 
      (b_groupVirgin8 + b_groupVirgin1)/2 - 
      (b_groupPregnant1 + b_groupPregnant8 + b_groupNone0) / 3
)
plot_post(flies3_post$contrast, xlab = "Virgin vs non-virgin groups")

La expresión

\[ \frac{\mu_3 + \mu_4}{2} - \frac{\mu_0 + \mu_1 + \mu_2}{3} = -\frac13 \mu_0 -\frac13 \mu_1 -\frac13 \mu_2 + \frac12 \mu_3 + \frac12 \mu_4 \]

es un ejemplo de un contraste. Un contraste es simplemente una combinación lineal de las medias de los grupos tal que la suma de los coeficientes sea 0. Muchas relaciones interesantes pueden investigarse usando contrastes, y el paquete brms incluye la función hypothesis() para ayudarnos a hacer esto. (Nota: debido a que incluimos sample_prior = TRUE en la llamada a brm() para este modelo, el gráfico a continuación muestra tanto las distribuciones previas como las posteriores para el contraste.)

h <-
  hypothesis(
    flies3_brm,
    "(groupVirgin8 + groupVirgin1) / 2 < 
      (groupPregnant1 + groupPregnant8 + groupNone0) / 3"
  )
h
plot(h)

8.1.4.4 Múltiples hipótesis a la vez

Incluso podemos probar múltiples hipótesis a la vez.

h2 <-
  hypothesis(
    flies3_brm,
    c("groupVirgin1 < (groupPregnant1 + groupPregnant8 + groupNone0) / 3",
      "groupVirgin8 < (groupPregnant1 + groupPregnant8 + groupNone0) / 3")
  )
h2
plot(h2)

8.1.5 Factores de Bayes

Supongamos que disponemos de datos \(D\) y queremos comparar una colección de modelos numerados por \(m=1,2,\dots\). Cada modelo \(m\) tiene sus propios parámetros \(\theta_m\), una verosimilitud

\[ p_m(D\mid \theta_m, m) \]

y una prior sobre esos parámetros

\[ p_m(\theta_m\mid m)\,. \]

Además asignamos a cada modelo una probabilidad a priori \(p(m)\).

Si consideramos el espacio conjunto de todos los \(\theta_m\) más el índice \(m\), la regla de Bayes en ese espacio se escribe

\[ p(\theta_1,\theta_2,\dots,m\mid D) = \frac{p(D\mid \theta_1,\theta_2,\dots,m)\;p(\theta_1,\theta_2,\dots,m)} {\displaystyle\sum_{m}\int p(D\mid \theta_1,\dots,m)\;p(\theta_1,\dots,m)\,d\theta_m}\,. \]

Bajo la suposición (jerárquica) de que, dado \(m\), los demás parámetros son independientes y sólo influyen a través de sus submodelos, esto se factoriza en

\[ p(\theta_1,\theta_2,\dots,m\mid D) = \frac{\prod_{m}\bigl[p_m(D\mid \theta_m,m)\,p_m(\theta_m\mid m)\,p(m)\bigr]} {\displaystyle\sum_{m}\int\prod_{m}\bigl[p_m(D\mid \theta_m,m)\,p_m(\theta_m\mid m)\,p(m)\bigr]\,d\theta_m}\,. \]

Sin embargo, si sólo estamos interesados en la probabilidad posterior del índice \(m\), marginalizamos todos los \(\theta_m\):

\[ p(m\mid D) = \frac{p(D\mid m)\,p(m)} {\sum_{m}p(D\mid m)\,p(m)}\,, \quad\text{donde}\quad p(D\mid m) = \int p_m(D\mid \theta_m,m)\;p_m(\theta_m\mid m)\,d\theta_m\,. \]

El factor de Bayes (BF) entre dos modelos \(m=1\) y \(m=2\) se define como la razón de sus evidencias:

\[ \mathrm{BF}_{1,2} = \frac{p(D\mid m=1)}{p(D\mid m=2)}\,. \]

Entonces las razones a posteriori y a priori de los modelos satisfacen

\[ \frac{p(m=1\mid D)}{p(m=2\mid D)} = \underbrace{\frac{p(m=1)}{p(m=2)}}_{\text{razón a priori}} \;\times\; \underbrace{\frac{p(D\mid m=1)}{p(D\mid m=2)}}_{\mathrm{BF}_{1,2}} \,. \]

Así, el factor de Bayes nos dice en cuánto multiplican las razones a priori al incorporar los datos. Una convención habitual es considerar evidencia “sustancial” a favor de \(m=1\) si \(\mathrm{BF}_{1,2}>3\), y a favor de \(m=2\) si \(\mathrm{BF}_{1,2}<\tfrac13\).

8.1.5.1 ¿Igualdad o desigualdad?

En el ejemplo anterior, expresamos nuestro contraste como una desigualdad. También podemos expresarlo como una igualdad. El resultado que obtenemos de hypothesis() es un poco diferente si lo hacemos así.

h3 <-
  hypothesis(
    flies3_brm,
    c("groupVirgin1 = (groupPregnant1 + groupPregnant8 + groupNone0) / 3",
      "groupVirgin8 = (groupPregnant1 + groupPregnant8 + groupNone0) / 3")
  )
h3
plot(h3)

• El IC es bidireccional en lugar de unidireccional.
• La Ratio de Evidencia se define de manera diferente.
  • Para las desigualdades, esta es la razón de las probabilidades posteriores de que la desigualdad sea verdadera vs. falsa.
  • Para las igualdades, esta es la razón de la densidad posterior (de la igualdad sostenida) a la densidad previa (Factor de Bayes). (Esto solo funciona si sample_prior = TRUE, ya que se requieren muestras previas para hacer el cálculo.)

8.2 ANOVA de varias vías

8.2.1 Rendimiento de Cultivo según el Método de Labranza y Fertilizante

Los datos en CalvinBayes::SplitPlotAgri provienen de un estudio agrícola en el que se utilizaron diferentes métodos de labranza y diferentes fertilizantes, y posteriormente se midió el rendimiento del cultivo (en bushels por acre) o (fanegas por acre).

gf_point(Yield ~ Fert | ~ Till, data = SplitPlotAgri, alpha = 0.4, size = 4)

Ajustamos dos modelos: sin/con interacción entre fertilizante y labranza. Estamos interesados en la pregunta: ¿Cómo usarías cada modelo para estimar el rendimiento medio al usar labranza de contorno (ridge) y fertilizante profundo (deep)?

fert1_brm <-
  brm(Yield ~ Till + Fert, data = SplitPlotAgri)

fert2_brm <-     
  brm(Yield ~ Till * Fert, data = SplitPlotAgri)

fert1_brm  

Nota que este modelo implica que la diferencia en rendimiento entre el uso de dos fertilizantes es la misma para cada uno de los tres métodos de labranza y la diferencia debida a los métodos de labranza es la misma para cada uno de los tres fertilizantes. Esto puede no ser una suposición razonable. Tal vez algunos fertilizantes funcionen mejor con ciertos métodos de labranza que con otros. El modelo 2 permite esto.

fert2_brm  

Como antes, podemos optar por ajustar el modelo sin una intersección. Esto produce una diferente parametrización del mismo modelo.

fert2a_brm <-   
  brm(Yield ~ 0 + Till * Fert, data = SplitPlotAgri)

fert2a_brm

8.2.2 Diseño de Parcela Dividida

El estudio utilizó 33 campos diferentes. Cada campo se dividió en 3 secciones y se aplicó un fertilizante diferente a cada una de las tres secciones. (Qué fertilizante se utilizó en qué sección se determinó al azar). Esto se llama un “diseño de parcela dividida” (incluso si se aplica a cosas que no son campos de cultivo).

Hubiera sido posible dividir cada campo en 9 subparcelas y usar todas las combinaciones de labranza y fertilizante, pero ese no fue el enfoque de este estudio. El método de labranza fue el mismo para todo el campo, probablemente porque era mucho más eficiente arar los campos de esta manera.

El gráfico a continuación indica que diferentes campos parecen tener rendimientos base diferentes, ya que los puntos asociados con un campo tienden a estar cerca de la parte superior o inferior de cada uno de los grupos de fertilizantes. Podemos agregar una variable adicional a nuestro modelo para manejar esta situación.

gf_point(Yield ~ Fert | ~ Till, data = SplitPlotAgri, alpha = 0.4, size = 4) %>%
  gf_line(group = ~Field)

fert3_brm <-
  # el uso de factor() es importante aquí porque los ids de campo son números
  # factor convierte esto en un factor (es decir, una variable nominal)
  brm(Yield ~ Till * Fert + factor(Field), data = SplitPlotAgri)

fert3_brm

Afortunadamente, en realidad no queremos este modelo. Ahora tenemos un ajuste para cada campo, y hubo 33 campos. Pero realmente no estamos interesados en predecir el rendimiento para un campo dado. Nuestro interés principal es en qué fertilizantes y métodos de labranza funcionan bien.

Si pensamos que la calidad del campo podría describirse mediante una distribución normal (u otra distribución), podríamos estar más interesados en los parámetros de esa distribución que en las estimaciones específicas para los campos particulares en este estudio. El tipo de modelo que queremos para esto se llama un modelo jerárquico o multinivel, y brm() facilita la descripción de dicho modelo.

Aquí hay una manera de pensar sobre tal modelo

• Cada campo tiene una productividad base.
• Las productividades base son normales con alguna media y desviación estándar que nos dicen sobre la distribución de la productividad entre campos. Nuestros 33 campos deberían ayudarnos a estimar esta distribución.
• Esa productividad base puede ajustarse hacia arriba o hacia abajo dependiendo del método de labranza y fertilizante utilizado.

En la jerga de brm(), el efecto del campo es ajustar la intersección, así que podemos escribirlo así:

fert4_brm <-
  brm(Yield ~ Till * Fert + (1 | Field), data = SplitPlotAgri)

Podemos ver en la salida a continuación que la variabilidad de parcela a parcela se estima con una desviación estándar de aproximadamente 8 a 15. Las estimaciones individuales de los campos están ocultas en este informe, pero puedes verlas si escribes stanfit(fert_brm).

fert4_brm

8.2.3 Medidas de precisión predictiva

Una forma de evaluar un modelo es mediante la precisión de sus predicciones. A veces nos interesa esta precisión en sí misma (por ejemplo, al evaluar pronósticos) y en otros casos la usamos para comparar varios modelos. Primero repasamos distintas maneras de cuantificar el error de predicción, y luego veremos cómo estimar esa precisión a partir de los datos.

8.2.3.1 Predicción puntual y funciones de puntuación

En la predicción puntual se estima un único valor \(\hat y_i\) para cada observación futura. Una función de puntuación común es el error cuadrático medio:

\[ \mathrm{MSE} = \frac1n\sum_{i=1}^n\bigl(y_i - E(y_i\mid\theta)\bigr)^2, \]

o su versión normalizada, dividiendo cada término por \(\text{Var}(y_i\mid\theta)\). Estas medidas son fáciles de calcular e interpretar, pero pierden adecuación si el modelo se aleja de la normalidad.

8.2.3.2 Predicción probabilística y reglas de puntuación

En la predicción probabilística se estudia la densidad predictiva completa \(p(\tilde y\mid\theta)\), y las reglas de puntuación (scoring rules) miden su calidad. Una buena regla es propia (incentiva reportar honestamente la creencia) y local (sanciona más ciertos errores). La única regla local y propia, salvo transformaciones afines, es la puntuación logarítmica:

\[ S(\tilde y;\theta) = \log p(\tilde y\mid\theta). \]

Ésta coincide, salvo constante, con el logverosimilitud cuando la varianza es constante.

8.2.3.3 Precisión predictiva para un solo punto

Supongamos que los datos verdaderos vienen de una distribución \(f(y)\), y que hemos obtenido la posterior \(p_{\rm post}(\theta)=p(\theta\mid y)\). La densidad predictiva para una futura observación \(\tilde y_i\) es

\[ p_{\rm post}(\tilde y_i) =\int p(\tilde y_i\mid\theta)\,p_{\rm post}(\theta)\,d\theta, \]

y su puntuación logarítmica es

\[ \log p_{\rm post}(\tilde y_i). \]

8.2.3.4 Promediado sobre futuros datos

Como \(\tilde y_i\) es desconocida, definimos la densidad predictiva esperada para un nuevo dato:

\[ \mathrm{elpd} \;=\;E_{f}\bigl[\log p_{\rm post}(\tilde y_i)\bigr] =\int \log p_{\rm post}(\tilde y_i)\,f(\tilde y_i)\,d\tilde y_i. \]

Para todo un conjunto de \(n\) puntos, usamos la suma punto a punto,

\[ \mathrm{elppd} =\sum_{i=1}^nE_{f}\bigl[\log p_{\rm post}(\tilde y_i)\bigr]. \tag{7.2} \]

También se puede calcular el valor esperado usando una estimación puntual \(\hat\theta(y)\):

\[ E_f\bigl[\log p(\tilde y\mid\hat\theta)\bigr]. \tag{7.3} \]

8.2.3.5 Evaluación práctica de la precisión predictiva

En la práctica no conocemos \(\theta\), pero disponemos de simulaciones posteriores \({\theta^s}_{s=1}^S\). Entonces la densidad predictiva punto a punto se estima como

\[ \mathrm{lppd} =\sum_{i=1}^n\log\int p(y_i\mid\theta)\,p_{\rm post}(\theta)\,d\theta \approx \sum_{i=1}^n\log\Bigl(\frac1S\sum_{s=1}^S p(y_i\mid\theta^s)\Bigr). \]

Este lppd sobreestima la elppd futura, por lo que luego aplicamos una corrección de sesgo.

8.2.4 Criterios de información y validación cruzada

Para comparar modelos según su precisión predictiva suele utilizarse la devianza, definido como $-2p(y),, $ donde \(\hat\theta\) es una estimación puntual ajustada al conjunto de datos \(y\).

Aunque un modelo ajustado rinde mejor sobre los datos de entrenamiento, su capacidad predictiva “out-of-sample” será en promedio menor. Para comparar modelos—incluso con distintas complejidades—necesitamos medidas que penalicen la tendencia de los modelos más grandes a sobreajustar.

8.2.4.1 Estimaciones con datos disponibles

1. Precisión dentro de muestra: usar el log-predictivo punto a punto (lppd) calculado anteriormente. Es rápido, pero sobreestima la precisión futura.

2. Precisión ajustada dentro de muestra: añadir una corrección al lppd para el número de parámetros (o “número efectivo” de parámetros). Esto da lugar a criterios como AIC, DIC y WAIC, que son aproximadamente insesgados en esperanza, pero no garantizan exactitud en cada caso.

3. Validación cruzada: dividir los datos en entrenamiento y prueba múltiples veces y promediar la precisión en los conjuntos de prueba. Evita el sobreajuste “in-sample” pero puede ser muy costosa computacionalmente; el LOO-CV (dejar uno fuera) implica en principio \(n\) ajustes del modelo, salvo atajos computacionales.

8.2.4.2 Akaike Information Criterion (AIC)

Cuando resumimos la inferencia de \(\theta\) con un estimador puntual \(\hat\theta\) (normalmente el MLE), la precisión predictiva out‐of‐sample se define como

\[ \mathrm{elpd}_{\hat\theta} =E_f\bigl[\log p(\tilde y\mid\hat\theta)\bigr]. \]

Como no podemos calcularlo directamente, usamos la log‐verosimilitud del ajuste a los datos observados y restamos un término de corrección por sobreajuste. Si \(k\) es el número de parámetros estimados,

\[ \widehat{\mathrm{elpd}}_{\rm AIC} =\log p\bigl(y\mid\hat\theta_{\rm mle}\bigr)\;-\;k, \]

y el AIC es

\[ \mathrm{AIC} =-2\,\log p\bigl(y\mid\hat\theta_{\rm mle}\bigr)\;+\;2k. \]

8.2.5 Deviance Information Criterion (DIC)

DIC adapta AIC al contexto bayesiano, reemplazando el MLE por la media posterior \(\hat\theta_{\rm Bayes}=E(\theta\mid y)\) y \(k\) por un número efectivo de parámetros \(p_{\rm DIC}\), definido como

\[ p_{\rm DIC} =2\Bigl[\log p\bigl(y\mid\hat\theta_{\rm Bayes}\bigr) -E_{\rm post}\bigl[\log p(y\mid\theta)\bigr]\Bigr]. \]

La corrección de sesgo da

\[ \widehat{\mathrm{elpd}}_{\rm DIC} =\log p\bigl(y\mid\hat\theta_{\rm Bayes}\bigr)\;-\;p_{\rm DIC}, \]

y el DIC se expresa como

\[ \mathrm{DIC} =-2\,\log p\bigl(y\mid\hat\theta_{\rm Bayes}\bigr)\;+\;2\,p_{\rm DIC}. \]

Una variante alternativa para \(p_{\rm DIC}\) es

\[ p_{\rm DIC}^{\rm alt} =2\,\text{Var}_{\rm post}\!\bigl[\log p(y\mid\theta)\bigr], \]

siendo siempre positiva.

8.2.6 Watanabe–Akaike Information Criterion (WAIC)

WAIC evalúa directamente la log‐densidad predictiva punto a punto,

\[ \mathrm{lppd} =\sum_{i=1}^n\log\Bigl(\tfrac1S\sum_{s=1}^S p(y_i\mid\theta^s)\Bigr), \]

y corrige el sobreajuste mediante un número efectivo de parámetros. Dos propuestas son:

\[ p_{\rm WAIC1} =2\sum_{i=1}^n\Bigl[\log E_{\rm post}p(y_i\mid\theta) -E_{\rm post}\log p(y_i\mid\theta)\Bigr], \]

\[ p_{\rm WAIC2} =\sum_{i=1}^n\text{Var}_{\rm post}\!\bigl[\log p(y_i\mid\theta)\bigr]. \]

Entonces

\[ \widehat{\mathrm{elpd}}_{\rm WAIC} =\mathrm{lppd}-p_{\rm WAIC}, \]

y el WAIC es \(-2\) veces este valor, quedando en la misma escala de la devianza. En la práctica \(p_{\rm WAIC2}\) suele ser más estable y coincidir mejor con LOO-CV (se ve más adelante).

8.2.7 Bayesian Information Criterion (BIC)

BIC penaliza más fuertemente según el tamaño de la muestra \(n\). Definido como

\[ \mathrm{BIC} =-2\,\log p\bigl(y\mid\hat\theta_{\rm mle}\bigr)\;+\;k\,\log n, \]

8.2.8 Validación cruzada “dejar-uno-fuera” (LOO-CV)

La validación cruzada bayesiana consiste en particionar los datos en un conjunto de entrenamiento \(y_{\rm train}\) y otro de prueba \(y_{\rm holdout}\), ajustar el modelo a \(y_{\rm train}\) para obtener la posterior

\[ p_{\rm train}(\theta)=p(\theta\mid y_{\rm train}), \]

y luego evaluar la densidad predictiva de los datos de prueba,

\[ \log p_{\rm train}(y_{\rm holdout}) =\log\int p(y_{\rm holdout}\mid\theta)\,p_{\rm train}(\theta)\,d\theta, \]

que se aproxima con \(S\) muestras \({\theta^s}\) de la posterior:

\[ \log\Bigl(\tfrac1S\sum_{s=1}^S p(y_{\rm holdout}\mid\theta^s)\Bigr). \]

En leave-one-out (LOO-CV), el conjunto de prueba contiene un solo dato \(y_i\), y repetimos este procedimiento \(n\) veces (cada vez dejando fuera otro dato), obteniendo \(n\) ajustes con muestras \(\theta^{(i)s}\). La log-densidad predictiva conjunta estimada es

\[ \mathrm{lppd}_{\rm loo\text{-}cv} =\sum_{i=1}^n \log p_{\rm post(-i)}(y_i) =\sum_{i=1}^n \log\Bigl(\tfrac1S\sum_{s=1}^S p(y_i\mid\theta^{(i)s})\Bigr). \]

Como cada predicción se condiciona en \(n-1\) datos, esta lppd tiende a subestimar el ajuste real para \(n\) datos completos. Para corregir el sesgo, se define

\[ b=\mathrm{lppd}-\overline{\mathrm{lppd}}_{-i}, \quad \overline{\mathrm{lppd}}_{-i} =\frac1n\sum_{i=1}^n\sum_{j=1}^n \log\Bigl(\tfrac1S\sum_{s=1}^S p(y_j\mid\theta^{(i)s})\Bigr), \]

y la versión corregida es

\[ \mathrm{lppd}_{\rm cloo\text{-}cv} =\mathrm{lppd}_{\rm loo\text{-}cv}+b. \]

Para comparar con otros criterios también se define un número efectivo de parámetros por

\[ p_{\rm loo\text{-}cv} =\mathrm{lppd}-\mathrm{lppd}_{\rm loo\text{-}cv}, \]

o usando la versión corregida.

Ventajas y limitaciones:

• Captura ajuste out-of-sample real, pero requiere reentrenar \(n\) veces (puede reducirse usando atajos).
• Asintóticamente (para \(n\to\infty\)) coincide con AIC, DIC y WAIC, pero para muestras pequeñas o modelos jerárquicos puede diferir.
• Requiere particionar datos en bloques condicionalmente independientes, lo cual puede ser problemático en series temporales, datos espaciales, etc.

8.2.9 Uso de loo

El paquete loo proporciona funciones para calcular estimaciones de WAIC y LOO de elpd (y sus contrapartes de criterios de información).

Aunque las definiciones son algo complejas, comparar modelos usando WAIC o LOO es relativamente sencillo. Según los autores del paquete loo, WAIC puede ser más rápido, pero LOO ofrece un mejor rendimiento.

library(loo)

waic(fert4_brm)

loo(fert4_brm)

A veces, el método de aproximación LOO-PSIS (importancia de muestreo suavizado de Pareto) no funciona bien y loo() recomienda reajustar algunos modelos desde cero. Esto se basa en el parámetro de forma (k) de la distribución de Pareto utilizada para suavizar las colas de la distribución posterior.

Permitamos que loo() vuelva a ajustar los modelos que considere necesario.

fert4_loo <- loo(fert4_brm, reloo = TRUE) # reajuste según sea necesario

En este caso, no hubo cambios significativos al reajustar los ocho modelos “problemáticos”.

fert4_loo
plot(fert4_loo)
fert4a_loo <- loo(fert4_brm)
plot(fert4a_loo)

Si tenemos múltiples modelos, podemos usar loo::compare() para compararlos según WAIC o LOO. Antes de hacerlo, agreguemos otro modelo a nuestra lista.

fert5_brm <- brm(Yield ~ Till + Fert + (1 | Field), data = SplitPlotAgri)

fert1_loo <- loo(fert1_brm)
fert2_loo <- loo(fert2_brm)
fert5_loo <- loo(fert5_brm)

Ahora podemos comparar nuestros cuatro modelos utilizando LOO:

compare(fert1_loo, fert2_loo, fert4_loo, fert5_loo)

Aspectos importantes a recordar:

1. Los elpd estimados y los criterios de información no son significativos por sí solos; solo son útiles para comparaciones.

2. Las comparaciones solo pueden hacerse entre modelos ajustados con los mismos datos, ya que los valores calculados dependen tanto del modelo como de los datos.

3. Todos estos métodos son aproximados. loo() y waic() proporcionan errores estándar además de las estimaciones.

4. p_loo (número efectivo de parámetros) también es una medida interesante. Si esta estimación no corresponde aproximadamente al número de parámetros libres en su modelo, suele ser señal de algún problema (posiblemente el modelo está mal especificado).