Criterios de Información

Introducción

Las medidas de precisión predictiva, conocidas como criterios de información, se definen en función de la devianza (deviance), que toma la forma:
\[ -2\,\log p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}})\,. \]

El factor \(-2\) se incluye para alinear esta medida con el estadístico utilizado en la prueba de razón de verosimilitudes (Likelihood Ratio Test), definido como:
\[ \lambda = -2\log\frac{\text{sup}_{\theta\in\Theta_0} L(\theta)}{\text{sup}_{\theta\in\Theta} L(\theta)}, \] donde \(L(\cdot)\) es la función de verosimilitud.

Sin embargo, seleccionar el modelo basado únicamente en la devianza más baja no es lo más adecuado. Es crucial realizar una corrección por la complejidad del modelo, considerando el número de parámetros estimados para evitar un ajuste excesivo (overfitting) y garantizar un equilibrio entre el ajuste y la capacidad predictiva del modelo.

Criterios de información

Existen varios métodos para estimar la precisión predictiva sin necesidad de usar datos fuera de la muestra. Entre ellos destacan los siguientes criterios de información:

Criterio de Información de la Devianza (DIC)

El DIC es una extensión Bayesiana del Criterio de Información de Akaike (AIC), el cual está definido como:
\[ \text{AIC} = -2\,\log p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{MLE}}) + 2k, \]
donde \(k\) es el número de parámetros del modelo.

En su versión Bayesiana, el DIC se define como:
\[ \text{DIC} = -2\,\log p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{Bayes}}) + 2p_{\text{DIC}}, \]
donde:

• \(\hat{\boldsymbol{\theta}}_{\text{Bayes}} = \mathsf{E}(\boldsymbol{\theta}\mid\boldsymbol{y}) \approx \frac{1}{B}\sum_{b=1}^B \boldsymbol{\theta}^{(b)}\) es la media posterior de \(\boldsymbol{\theta}\).
• \(p_{\text{DIC}}\) es el número efectivo de parámetros, dado por: \[ p_{\text{DIC}} = 2\left(\log p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{Bayes}}) - \mathsf{E}(\log p(\boldsymbol{y}\mid\boldsymbol{\theta}) \mid \boldsymbol{y})\right), \] lo cual se aproxima como: \[ p_{\text{DIC}} \approx 2\left(\log p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{Bayes}}) - \frac{1}{B}\sum_{b=1}^B \log p(\boldsymbol{y}\mid\boldsymbol{\theta}^{(b)})\right). \]

Referencias:

• Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4), 583-639.

Criterio de Información de Watanabe-Akaike (WAIC)

El WAIC es un criterio completamente Bayesiano para medir la precisión predictiva y se define como:
\[ \text{WAIC} = -2\,\text{lppd} + 2p_{\text{WAIC}}, \]
donde:

• \(\text{lppd}\) es la densidad predictiva puntual logarítmica: \[ \text{lppd} = \sum_{i=1}^n \log \int_\Theta p(y_i \mid \boldsymbol{\theta}) p(\boldsymbol{\theta}\mid\boldsymbol{y}) \, \mathrm{d}\boldsymbol{\theta} \approx \sum_{i=1}^n \log\left(\frac{1}{B}\sum_{b=1}^B p(y_i \mid \boldsymbol{\theta}^{(b)})\right). \]
• \(p_{\text{WAIC}}\) es el número efectivo de parámetros, dado por: \[ p_{\text{WAIC}} = 2\sum_{i=1}^n\left(\log \mathsf{E}(p(y_i\mid\boldsymbol{\theta}) \mid \boldsymbol{y}) - \mathsf{E}(\log p(y_i\mid\boldsymbol{\theta}) \mid \boldsymbol{y})\right), \] lo cual se aproxima como: \[ p_{\text{WAIC}} \approx 2\sum_{i=1}^n\left(\log\left(\frac{1}{B}\sum_{b=1}^B p(y_i \mid \boldsymbol{\theta}^{(b)})\right) - \frac{1}{B}\sum_{b=1}^B \log p(y_i \mid \boldsymbol{\theta}^{(b)})\right). \]

Referencias:

• Watanabe, S., & Opper, M. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11(12).
• Gelman, A., Hwang, J., & Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24(6), 997-1016.

Criterio de Información Bayesiano (BIC)

El BIC penaliza la complejidad del modelo en función del tamaño de la muestra, y se define como: \[ \text{BIC} = -2\,\log p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{Bayes}}) + k\,\log(n), \]
donde: - \(k\) es el número de parámetros del modelo. - \(n\) es el tamaño de la muestra.

Referencias:

• Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464.

Ejemplo: Puntajes de Matemáticas

La base de datos contiene información proveniente de una muestra aleatoria simple de los estudiantes que presentaron la Prueba Saber 11 en el segundo semestre de 2023.

La prueba de matemáticas está construida sobre una escala de 0 a 100 (sin decimales), con un puntaje promedio de 50 y una desviación estándar de 10 puntos.

El objetivo es construir un modelo para el puntaje de matemáticas a nivel nacional, utilizando como datos de entrenamiento los resultados del segundo semestre de 2023, con el propósito de realizar inferencias sobre la población estudiantil tanto a nivel nacional como departamental.

Por esta razón, se considera el departamento de residencia del estudiante como variable de agrupamiento.

Los datos son de acceso público y pueden consultarse en el siguiente enlace.

Estructura de los datos

• \(y_{i,j}\): puntaje de matemáticas del estudiante \(i\) y en departamento \(j\).
• \(n_j\,\,\): número de estudiantes en el departamento \(j\).
• \(\bar{y}_j\,\,\): promedio muestral del departamento \(j\).
• \(s^2_j\,\,\): varianza muestral del departamento \(j\).

Tratamiento de datos

# Datos
dat <- read.csv("SB11_20232_muestra.txt", sep=";")
dat <- dat[order(dat$ESTU_COD_RESIDE_MCPIO), ]

# Dimensión de la base
dim(dat)

## [1] 5511   83

# Distribución de frecuencias
table(dat$ESTU_DEPTO_RESIDE)

## 
##        AMAZONAS       ANTIOQUIA          ARAUCA       ATLANTICO          BOGOTÁ 
##               8             758              34             342             776 
##         BOLIVAR          BOYACA          CALDAS         CAQUETA        CASANARE 
##             275             173              78              42              53 
##           CAUCA           CESAR           CHOCO         CORDOBA    CUNDINAMARCA 
##             126             149              55             218             389 
##         GUAINIA        GUAVIARE           HUILA      LA GUAJIRA       MAGDALENA 
##               4              15             126             100             200 
##            META          NARIÑO NORTE SANTANDER        PUTUMAYO         QUINDIO 
##             122             145             147              38              45 
##       RISARALDA      SAN ANDRES       SANTANDER           SUCRE          TOLIMA 
##             125               9             274             107             170 
##           VALLE          VAUPES         VICHADA 
##             401               3               4

# paquetes
suppressMessages(suppressWarnings(library(dplyr))) 
suppressMessages(suppressWarnings(library(ggplot2))) 

# m : número de grupos (departamentos)
(m <- length(table(dat$ESTU_DEPTO_RESIDE)))

## [1] 33

# n : número de individuos (estudiantes)
(n <- sum(table(dat$ESTU_DEPTO_RESIDE)))

## [1] 5511

# Tratamiento de datos
# y  : puntaje de los estudiantes (c)
# Y  : puntaje de los estudiantes (list)
# g  : identificador secuencial de los departamentos (c)
# nj : número de estudiantes por departamento (c)
# yb : promedios por departamento (c)
# s2 : varianzas por departamento (c)

y <- dat$PUNT_MATEMATICAS
Y <- vector(mode = "list", length = m)
g <- rep(NA, n)
for (j in 1:m) {
  idx <- dat$ESTU_COD_RESIDE_DEPTO == sort(unique(dat$ESTU_COD_RESIDE_DEPTO))[j]
  g[idx] <- j
  Y[[j]] <- y[idx]
}

# Tabla
estadisticos <- dat %>% 
  group_by(ESTU_COD_RESIDE_DEPTO) %>% 
  summarise(
    codigo = first(ESTU_COD_RESIDE_DEPTO),
    nombre = first(ESTU_DEPTO_RESIDE),
    nj = n(), 
    yb = mean(PUNT_MATEMATICAS), 
    s2 = var(PUNT_MATEMATICAS)
  ) %>% 
  ungroup() %>% 
  select(-ESTU_COD_RESIDE_DEPTO)

# Encabezado de la tabla
head(estadisticos, n = 10)

## # A tibble: 10 × 5
##    codigo nombre       nj    yb    s2
##     <int> <chr>     <int> <dbl> <dbl>
##  1      5 ANTIOQUIA   758  49.9  140.
##  2      8 ATLANTICO   342  50.1  140.
##  3     11 BOGOTÁ      776  55.1  144.
##  4     13 BOLIVAR     275  46.8  177.
##  5     15 BOYACA      173  54.0  150.
##  6     17 CALDAS       78  52.2  119.
##  7     18 CAQUETA      42  50.1  167.
##  8     19 CAUCA       126  48.0  164.
##  9     20 CESAR       149  50.5  137.
## 10     23 CORDOBA     218  49.1  182.

# Tamaños de muestra
nj <- estadisticos$nj

# Estadísticos suficientes
yb <- estadisticos$yb
s2 <- estadisticos$s2

Modelamiento

• \(\textsf{M}_1\): Modelo jerárquico Normal con medias específicas.
• \(\textsf{M}_2\): Modelo jerárquico Normal con medias y varianzas específicas.

MCMC1 <- function(B, y, nj, yb, s2, mu0, g20, eta0, t20, nu0, s20) {
  # frecuencia para imprimir progreso
  ncat <- floor(B / 10) 
  
  # número total de observaciones y grupos
  n <- sum(nj)
  m <- length(nj)
  
  # almacenamiento de la cadena
  THETA <- matrix(NA, nrow = B, ncol = m + 4)
  
  # valores iniciales
  theta <- yb
  sig2  <- mean(s2)
  mu    <- mean(theta)
  tau2  <- var(theta)
  
  # cadena de MCMC
  for (b in 1:B) {
    
    # actualizar theta_j
    v_theta <- 1 / (1 / tau2 + nj / sig2)
    m_theta <- v_theta * (mu / tau2 + nj * yb / sig2)
    theta  <- rnorm(m, mean = m_theta, sd = sqrt(v_theta))
    
    # actualizar sigma^2
    a_sig2 <- 0.5 * (nu0 + n)
    b_sig2 <- 0.5 * (nu0 * s20 + sum((nj - 1) * s2 + nj * (yb - theta)^2))
    sig2 <- 1 / rgamma(1, shape = a_sig2, rate = b_sig2)
    
    # actualizar mu
    v_mu <- 1 / (1 / g20 + m / tau2)
    m_mu <- v_mu * (mu0 / g20 + m * mean(theta) / tau2)
    mu  <- rnorm(1, mean = m_mu, sd = sqrt(v_mu))
    
    # actualizar tau^2
    a_tau2 <- 0.5 * (eta0 + m)
    b_tau2 <- 0.5 * (eta0 * t20 + (m - 1) * var(theta) + m * (mean(theta) - mu)^2)
    tau2 <- 1 / rgamma(1, shape = a_tau2, rate = b_tau2)
    
    # log-verosimilitud
    ll <- sum(dnorm(x = y, mean = rep(theta, nj), sd = sqrt(sig2), log = TRUE))
    
    # almacenar iteración
    THETA[b, ] <- c(theta, sig2, mu, tau2, ll)
    
    # imprimir progreso
    if (b %% ncat == 0) cat(100 * round(b / B, 1), "% completado ... \n", sep = "")
  }
  
  # salida
  colnames(THETA) <- c(paste0("theta", 1:m), "sig2", "mu", "tau2", "ll")
  THETA <- as.data.frame(THETA)
  return(list(THETA = THETA))
}

MCMC2 <- function(B, y, nj, yb, s2, mu0, g20, eta0, t20, lam0, al0, be0, nus0) {
  # Ajustes
  ncat <- floor(B / 10)
  
  # Tamaños
  n <- sum(nj)
  m <- length(nj)
  
  # Almacenamiento
  THETA <- matrix(NA, nrow = B, ncol = 2 * m + 5)
  
  # Valores iniciales
  theta <- yb
  sig2  <- s2           # varianzas específicas por grupo
  mu    <- mean(theta)  
  tau2  <- var(theta)
  nu    <- 1
  ups2  <- 100.         # varianza global
  
  # Cadena MCMC
  for (b in 1:B) {
    
    # Actualizar theta_j
    v_theta <- 1 / (1 / tau2 + nj / sig2)
    m_theta <- v_theta * (mu / tau2 + nj * yb / sig2)
    theta   <- rnorm(m, mean = m_theta, sd = sqrt(v_theta))
    
    # Actualizar sigma_j^2
    a_sig2 <- 0.5 * (nu + nj)
    b_sig2 <- 0.5 * (nu * ups2 + (nj - 1) * s2 + nj * (yb - theta)^2)
    sig2 <- 1 / rgamma(m, shape = a_sig2, rate  = b_sig2)
    
    # Actualizar mu
    v_mu <- 1 / (1 / g20 + m / tau2)
    m_mu <- v_mu * (mu0 / g20 + m * mean(theta) / tau2)
    mu   <- rnorm(1, mean = m_mu, sd = sqrt(v_mu))
    
    # Actualizar tau^2
    a_tau2 <- 0.5 * (eta0 + m)
    b_tau2 <- 0.5 * (eta0 * t20 + (m - 1) * var(theta) + m * (mean(theta) - mu)^2)
    tau2 <- 1 / rgamma(1, shape = a_tau2, rate  = b_tau2)
    
    # Actualizar nu
    lpnu <- 0.5 * m * nus0 * log(0.5 * nus0 * ups2) -
            m * lgamma(0.5 * nus0) -
            0.5 * nus0 * sum(log(sig2)) -
            nus0 * (lam0 + 0.5 * ups2 * sum(1 / sig2))
    
    nu <- sample(x = nus0, size = 1, prob = exp(lpnu - max(lpnu)))
    
    # Actualizar ups2 (sigma^2 global)
    a_ups2 <- al0 + 0.5 * m * nu
    b_ups2 <- be0 + 0.5 * nu * sum(1 / sig2)
    ups2 <- rgamma(1, shape = a_ups2, rate = b_ups2)
    
    # Log-verosimilitud
    ll <- sum(dnorm(x = y, mean = rep(theta, nj), sd = sqrt(rep(sig2, nj)), log = TRUE))
    
    # Almacenar resultados
    THETA[b, ] <- c(theta, sig2, mu, tau2, nu, ups2, ll)
    
    # Progreso
    if (b %% ncat == 0) {
      cat(100 * round(b / B, 1), "% completado ... \n", sep = "")
    }
  }
  
  # Salida final
  colnames(THETA) <- c(paste0("theta", 1:m), paste0("sig2", 1:m), "mu", "tau2", "nu", "ups2", "ll")
  THETA <- as.data.frame(THETA)
  return(list(THETA = THETA))
}

Ajuste del modelo 1:

# Hiperparámetros
mu0  <- 50 
g20  <- 10^2
eta0 <- 1  
t20  <- 10^2
nu0  <- 1  
s20  <- 10^2

# ajuste del modelo
set.seed(123)
chain1 <- MCMC1(B = 10000, y, nj, yb, s2, mu0, g20, eta0, t20, nu0, s20)

## 10% completado ... 
## 20% completado ... 
## 30% completado ... 
## 40% completado ... 
## 50% completado ... 
## 60% completado ... 
## 70% completado ... 
## 80% completado ... 
## 90% completado ... 
## 100% completado ...

Ajuste del modelo 2:

# Hiperparámetros
mu0  <- 50 
g20  <- 10^2
eta0 <- 1  
t20  <- 10^2
lam0 <- 1  
al0  <- 1
be0  <- 1/10^2 
nus0 <- 1:50  # rango para p(nu | rest)

# ajuste del modelo
set.seed(123)
chain2 <- MCMC2(B = 10000, y, nj, yb, s2, mu0, g20, eta0, t20, lam0, al0, be0, nus0)

## 10% completado ... 
## 20% completado ... 
## 30% completado ... 
## 40% completado ... 
## 50% completado ... 
## 60% completado ... 
## 70% completado ... 
## 80% completado ... 
## 90% completado ... 
## 100% completado ...

Log-verosimilitud

Criterios de información

# DIC: Modelo 1
LL1        <- chain1$THETA$ll

theta_hat  <- colMeans(chain1$THETA[, 1:m])
sigma2_hat <- mean(chain1$THETA$sig2)

lpy_m1     <- sum(dnorm(x = y, mean = rep(theta_hat, nj), sd = sqrt(sigma2_hat), log = TRUE))

pDIC_m1    <-  2 * (lpy_m1 - mean(LL1))
dic_m1     <- -2 * lpy_m1 + 2 * pDIC_m1

# DIC: Modelo 2
LL2        <- chain2$THETA$ll

theta_hat  <- colMeans(chain2$THETA[,1:m])
sigma2_hat <- colMeans(chain2$THETA[,(m+1):(2*m)])

lpy_m2     <- sum(dnorm(x = y, mean = rep(theta_hat, nj), sd = sqrt(rep(sigma2_hat, nj)), log = T))

pDIC_m2    <-  2*(lpy_m2 - mean(LL2))
dic_m2     <- -2*lpy_m2 + 2*pDIC_m2

# WAIC
lppd_m1  <- 0
pWAIC_m1 <- 0
for (i in 1:n) {
  # lppd
  tmp1    <- dnorm(x = y[i], mean = chain1$THETA[, g[i]], sd = sqrt(chain1$THETA$sig2))
  lppd_m1 <- lppd_m1 + log(mean(tmp1))

  # pWAIC
  tmp2 <- dnorm(x = y[i], mean = chain1$THETA[, g[i]], sd = sqrt(chain1$THETA$sig2), log = TRUE)
  pWAIC_m1 <- pWAIC_m1 + 2 * (log(mean(tmp1)) - mean(tmp2))
}

waic_m1 <- -2 * lppd_m1 + 2 * pWAIC_m1

# WAIC
lppd_m2  <- 0
pWAIC_m2 <- 0
for (i in 1:n) {
  # lppd
  tmp1 <- dnorm(x = y[i], mean = chain2$THETA[,g[i]], sd = sqrt(chain2$THETA[,m+g[i]]))
  lppd_m2 <- lppd_m2 + log(mean(tmp1))
  # pWAIC
  tmp2 <- dnorm(x = y[i], mean = chain2$THETA[,g[i]], sd = sqrt(chain2$THETA[,m+g[i]]), log = T)
  pWAIC_m2 <- pWAIC_m2 + 2*(log(mean(tmp1)) - mean(tmp2))
}
waic_m2 <- -2*lppd_m2 + 2*pWAIC_m2

	Modelo 1	Modelo 2
lp	-21545.13	-21528.09
pDIC	28.50	52.52
DIC	43147.26	43161.22
lppd	-21545.34	-21530.15
pWAIC	28.08	48.39
WAIC	43146.84	43157.09

Ejercicios

• Considere el modelo \(y_i \mid \theta, \sigma^2 \overset{\text{ind}}{\sim} \textsf{t}_\kappa(\theta_j, \sigma^2)\), para \(i = 1, \dots, n\), donde \(\kappa \sim \textsf{U}\{1, 2, \dots, \nu_0\}\), con \(\nu_0\) como hiperparámetro. El resto del modelo se especifica de forma análoga al modelo normal jerárquico con medias específicas por grupo y varianza común.

  1. Reescriba la distribución muestral del modelo utilizando variables auxiliares, mediante la representación de la distribución t como una composición de una distribución Normal con una distribución Gamma Inversa (ver los ejercicios del “Muestreador de Gibbs”).
  2. Derive las distribuciones condicionales completas de todos los parámetros del modelo.
  3. Ajuste el modelo a los datos de la prueba Saber 11 del segundo semestre de 2023.
  4. Identifique los valores atípicos (outliers) utilizando las variables auxiliares.
  5. Realice toda la inferencia que se llevó a cabo en el caso del modelo normal jerárquico con medias específicas por grupo y varianza común.
  6. Interprete y compare los resultados obtenidos con los del modelo normal jerárquico.
  7. Calcule todos los criterios de información y compare los resultados con los modelos normales.

• Repetir el ejercicio anterior para el modelo \(y_i \mid \theta, \sigma^2 \overset{\text{ind}}{\sim} \textsf{t}_\kappa(\theta_j, \sigma^2_j)\), para \(i = 1, \dots, n\).

Referencias

Hoff, P. D. (2009). A First Course in Bayesian Statistical Methods. Springer New York.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). Chapman & Hall/CRC.