Introducción
El tipo más simple de datos multinivel contempla dos niveles: grupos y unidades dentro de los grupos.
Se denota con \(y_{i,j}\) la observación correspondiente a la unidad \(i\) del grupo \(j\), donde \(i = 1,\ldots,n_j\) y \(j = 1,\ldots,m\), siendo \(m\) el número total de grupos y \(n = \sum_{j=1}^m n_j\) el tamaño total de la muestra.
El conjunto de datos se representa como \(\boldsymbol{y} = (\boldsymbol{y}_1,\ldots,\boldsymbol{y}_m)\), donde \(\boldsymbol{y}_j = (y_{1,j},\ldots,y_{n_j,j})\) contiene las observaciones del grupo \(j\), para \(j = 1,\ldots,m\).
Modelo Normal multinivel
Un modelo ampliamente utilizado para describir la heterogeneidad de las medias entre distintas poblaciones es el modelo jerárquico Normal, en el cual la variabilidad dentro y entre los grupos se representa mediante distribuciones Normales:
Caracterización dentro de los grupos:
\[ y_{i,j} \mid \theta_j, \sigma^2 \overset{\text{ind}}{\sim} \textsf{N}(\theta_j, \sigma^2) \]
Caracterización entre los grupos:
\[ \theta_j \mid \mu, \tau^2 \overset{\text{iid}}{\sim} \textsf{N}(\mu, \tau^2) \]
Distribución previa conjunta:
\[ p(\sigma^2, \mu, \tau^2) = p(\sigma^2)\, p(\mu)\, p(\tau^2) \]
con especificaciones:
\[ \sigma^2 \sim \textsf{GI}\left(\tfrac{\nu_0}{2}, \tfrac{\nu_0\, \sigma^2_0}{2}\right), \qquad \mu \sim \textsf{N}(\mu_0, \gamma_0^2), \qquad \tau^2 \sim \textsf{GI}\left(\tfrac{\eta_0}{2}, \tfrac{\eta_0\, \tau^2_0}{2}\right) \]
Los parámetros del modelo son \(\boldsymbol{\theta} = (\theta_1, \ldots, \theta_m, \sigma^2, \mu, \tau^2)\).
Los hiperparámetros del modelo son \((\nu_0, \sigma^2_0, \mu_0, \gamma_0^2, \eta_0, \tau^2_0)\).
Entrenamiento
La estimación de los parámetros del modelo se puede hacer por medio de un muestreador de Gibbs para obtener muestras de la distribución posterior \(p(\boldsymbol{\theta}\mid\boldsymbol{y})\).
(Ejercicio.) Distribución posterior: \[ \begin{aligned} p(\boldsymbol{\theta} \mid \boldsymbol{y}) &\propto p(\boldsymbol{y} \mid \boldsymbol{\theta})\, p(\boldsymbol{\theta}) \\ &=\prod_{j=1}^m\prod_{i=1}^{n_j} \textsf{N}\left(y_{i, j} \mid \theta_{j}, \sigma^{2}\right) \\ &\quad\quad\times \prod_{j=1}^m \textsf{N}\left(\theta_{j} \mid \mu, \tau^{2}\right) \times \textsf{GI}\left(\sigma^{2} \mid \tfrac{\nu_{0}}{2}, \tfrac{\nu_{0}\,\sigma_{0}^{2}}{2} \right) \\ &\quad\quad\quad\quad\times \textsf{N}\left(\mu \mid \mu_{0}, \gamma_{0}^{2}\right) \times \textsf{GI}\left(\tau^{2} \mid \tfrac{\eta_{0}}{2}, \tfrac{\eta_{0}\,\tau_{0}^{2}}{2}\right) \end{aligned} \]
(Ejercicio.) Distribuciones condicionales completas: \[ \begin{aligned} \theta_{j} \mid - &\sim \textsf{N}\left( \frac{\mu / \tau^{2} + n_{j} \bar{y}_{j} / \sigma^{2}}{1 / \tau^{2} + n_{j} /\sigma^{2}}, \frac{1}{1 / \tau^{2} + n_{j} /\sigma^{2}}\right) \\ \sigma^{2} \mid - &\sim\textsf{GI}\left( \frac{\nu_{0}+\sum_{j=1}^m n_{j}}{2}, \frac{\nu_{0} \sigma_{0}^{2}+\sum_{j=1}^m \sum_{i=1}^{n_j}\left(y_{i, j}-\theta_{j}\right)^{2}}{2}\right) \\ \mu \mid - &\sim\textsf{N}\left( \frac{\mu_{0} / \gamma_{0}^{2} + m \bar{\theta} / \tau^{2}}{1 / \gamma_{0}^{2} + m / \tau^{2}}, \frac{1}{1 / \gamma_{0}^{2} + m / \tau^{2}}\right) \\ \tau^{2} \mid - &\sim\textsf{ GI }\left(\frac{\eta_{0}+m}{2}, \frac{\eta_{0} \tau_{0}^{2}+\sum_{j=1}^m\left(\theta_{j}-\mu\right)^{2}}{2}\right) \end{aligned} \]
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), ]
# Dimensiones
dim(dat)## [1] 5511 83##
## 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)))
suppressMessages(suppressWarnings(library(corrplot)))## [1] 33## [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)## # 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.Análisis exploratorio
# Cargar shapefile de departamentos de Colombia
# Fuente: https://sites.google.com/site/seriescol/shapes
suppressMessages(suppressWarnings(library(sf)))
suppressMessages(suppressWarnings(library(dplyr)))
suppressMessages(suppressWarnings(library(ggplot2)))
shp <- st_read("depto.shp", quiet = TRUE)
# Preparar datos de promedio por departamento
dat_map <- estadisticos[, c("codigo", "yb")]
pd <- dat_map
pd$codigo <- as.character(pd$codigo)
# Normalizar códigos con ceros a la izquierda
pd$codigo[pd$codigo == "5"] <- "05"
pd$codigo[pd$codigo == "8"] <- "08"
# Renombrar columnas
colnames(pd) <- c("DPTO", "Media")
# Unir shapefile con datos y graficar
shp %>%
inner_join(pd, by = "DPTO") %>%
select(DPTO, Media, geometry) %>%
ggplot() +
geom_sf(aes(fill = Media), size = 0.125, color = "#b2b2b2") +
theme_bw() +
labs(
title = "",
x = "Longitud",
y = "Latitud",
fill = "Media"
)
# Ranking basado en el promedio muestral
par(mfrow = c(1, 1), mar = c(4, 10, 1.5, 1), mgp = c(2.5, 0.75, 0))
# Inicializar gráfico vacío
plot(x = c(0, 100), y = c(1, m), type = "n",
xlab = "Puntaje", ylab = "",
main = expression(italic("Ranking (promedio muestral)")),
yaxt = "n")
# Líneas horizontales de referencia
abline(h = 1:m, col = "lightgray", lwd = 1)
# Línea vertical en el promedio nacional (50)
abline(v = 50, col = "gray", lwd = 3)
# Dibujar puntos individuales por departamento según el orden del promedio
for (l in 1:m) {
j <- order(yb)[l]
points(x = Y[[j]], y = rep(l, nj[j]), pch = 16, cex = 0.4)
}
# Añadir puntos del promedio muestral
lines(x = yb[order(yb)], y = 1:m, type = "p", col = 4, pch = 16, cex = 1.1)
# Conectar promedios con líneas
lines(x = yb[order(yb)], y = 1:m, type = "l", col = adjustcolor(4, alpha.f = 0.3))
# Etiquetas con nombres de los departamentos
axis(side = 2, at = 1:m, labels = estadisticos$nombre[order(yb)], las = 2)
# Configurar dos paneles lado a lado
par(mfrow = c(1, 2), mar = c(3, 3, 1.5, 1), mgp = c(1.75, 0.75, 0))
# Histograma del promedio por grupo
hist(yb, freq = FALSE,
nclass = 15,
main = "",
xlab = "Promedio",
ylab = "Densidad",
border = adjustcolor(4, alpha.f = 0),
col = adjustcolor(4, alpha.f = 0.3))
abline(v = mean(y), col = "gray", lwd = 3)
# Diagrama de dispersión: tamaño del grupo vs. promedio
plot(nj, yb,
xlab = "Tamaño del grupo",
ylab = "Promedio",
pch = 16, cex = 1.2,
col = adjustcolor(4, alpha.f = 0.6))
abline(h = mean(y), col = "gray", lwd = 3)
Es habitual que los grupos con promedios muestrales extremadamente altos o bajos correspondan a aquellos con tamaños muestrales reducidos.
Distribución previa
Teniendo en cuenta la información de la prueba, el modelo se ajusta utilizando los siguientes hiperparámetros: \[ \mu_0 = 50\,,\qquad \gamma^2_0 = 10^2\,,\qquad \eta_0 = 1\,,\qquad \tau^2_0 = 10^2\,,\qquad \nu_0 = 1\,,\qquad \sigma^2_0 = 10^2\,. \]
Ajuste del modelo
Inicializar \(\theta_1,\ldots,\theta_m\), \(\sigma^2\), \(\mu\), \(\tau^2\).
Repetir para \(b = 1, \ldots, B\):
- Actualizar \(\theta_j \, \sim p(\theta_j \mid -)\).
- Actualizar \(\sigma^2 \sim p(\sigma^2 \mid -)\).
- Actualizar \(\mu \,\,\sim p(\mu \mid -)\).
- Actualizar \(\tau^2 \sim p(\tau^2 \mid -)\).
- Calcular log-verosimilitud con \(\boldsymbol{y} \mid \theta_1,\ldots,\theta_m, \sigma^2\).
- Almacenar parámetros y log-verosimilitud.
Devolver la cadena simulada.
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))
}# 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 ...Convergencia
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5714 9555 10000 9649 10000 10745# Coeficiente de variación de Monte Carlo (%)
EEMC1 <- apply(X = chain1$THETA, MARGIN = 2, FUN = sd)/sqrt(neff1)
CVMC1 <- 100*abs(EEMC1/colMeans(chain1$THETA))
round(summary(CVMC1), 4)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0002 0.0168 0.0203 0.0383 0.0374 0.3989
Inferencia
Se define \(\eta\) como la proporción de intravarianza: \(\eta = 100\cdot\frac{\sigma^2}{\sigma^2+\tau^2}\%\).
- El 90.0% de la variabilidad total en los puntajes se explica por las
diferencias dentro de los departamentos.
- La media global de los puntajes es cercana a 50.
- El 99.7% de los puntajes individuales dentro de los departamentos se
encuentran a una distancia máxima de \(3
\times 12.1 = 36.3\) puntos.
- El 99.7% de los promedios de puntajes entre departamentos están separados por una distancia máxima de \(3 \times 4.0 = 12.0\) puntos.
| Media | CV(%) | Q2.5% | Q97.5% | |
|---|---|---|---|---|
| eta | 90.050 | 2.942 | 83.816 | 94.132 |
| mu | 49.865 | 1.554 | 48.331 | 51.357 |
| sig | 12.099 | 0.961 | 11.876 | 12.329 |
| tau | 3.997 | 14.629 | 3.019 | 5.322 |
Ranking
Contracción
La contracción (shrinkage) es un fenómeno donde las estimaciones de ciertos parámetros tienden a “acercarse” o “ajustarse” hacia un valor central o promedio. Este efecto ocurre generalmente debido a la influencia de la previa o a la estructura jerárquica del modelo.
La contracción ayuda a proporcionar estimaciones más moderadas y estables, especialmente en escenarios donde los datos son limitados.
Este ajuste mejora la robustez de los resultados al equilibrar la información observada con la información previa o el contexto global del modelo.
La contracción genera dependencia entre los grupos al asumir que sus parámetros comparten una distribución común. Esto provoca que las estimaciones de los grupos con poca información se ajusten hacia el promedio global, mientras que los grupos con mayor cantidad de información se apoyen más en sus propios datos.
Agrupamiento
Se aplica un algoritmo de segmentación (e.g., k-means) sobre las medias posteriores de \((\theta_1, \sigma^2_1), \ldots, (\theta_m, \sigma^2_m)\).
Para seleccionar el número óptimo de grupos, se implemente un k-means con distintos valores de \(K\) (e.g., 2 a 15) y se calcula la suma de cuadrados intra-cluster (WSS, within sum of squares). A partir de la reducción relativa \(\Delta_k = (\text{WSS}_{k-1} - \text{WSS}_{k})/\text{WSS}_{k-1}\), se elige el menor valor de \(k\) tal que \(\Delta_k < \varepsilon\) (e.g., \(\varepsilon = 0.2\)), es decir, el punto en que añadir más grupos no mejora sustancialmente el ajuste.
# Número de grupos (departamentos)
m <- length(table(dat$ESTU_DEPTO_RESIDE))
# Número de iteraciones y departamentos
B <- nrow(chain1$THETA[ , 1:m])
# Rango de valores de k
k_range <- 2:15
# umbral de mejora relativa
eps <- 0.2
# Media posterior de cada grupo
theta_hat <- colMeans(chain1$THETA[ , 1:m])
# Inicializar almacenamiento para WSS y clusters
wss <- numeric(length(k_range))
clusters_list <- vector("list", length(k_range))
# Evaluar k-means para cada valor de k en k_range
for (i in seq_along(k_range)) {
k <- k_range[i]
km <- kmeans(theta_hat, centers = k, nstart = 10)
wss[i] <- km$tot.withinss
clusters_list[[i]] <- km$cluster
}
# Calcular reducción relativa de WSS
delta <- abs(diff(wss)) / wss[-length(wss)]
best_k <- which(delta < eps)
# Seleccionar número óptimo de clusters
if (length(best_k) == 0) {
xi_hat <- clusters_list[[1]] # usar k = min(k_range) si no hay mejora significativa
} else {
xi_hat <- clusters_list[[min(best_k) + 1]]
}
# Número de clusters
(K <- length(table(xi_hat)))## [1] 8# Leer shapefile
shp <- st_read("depto.shp", quiet = TRUE)
# Preparar data frame con códigos y clusters
pd <- estadisticos[, c("codigo")]
pd$codigo <- as.character(pd$codigo)
# Normalizar códigos con ceros a la izquierda si es necesario
pd$codigo[pd$codigo == "5"] <- "05"
pd$codigo[pd$codigo == "8"] <- "08"
# Agregar clusters
pd$Cluster <- xi_hat
colnames(pd)[1] <- "DPTO" # renombrar columna
# Generar colores interpolados
colores <- c(
"#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00",
"#FFFF33", "#A65628", "#F781BF", "#999999", "#66C2A5",
"#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F",
"#E5C494", "#B3B3B3", "#1B9E77", "#D95F02", "#7570B3"
)[1:K]
# Unir shapefile con datos y graficar
shp %>%
inner_join(pd, by = "DPTO") %>%
select(DPTO, Cluster, geometry, nombre = NOMBRE_DPT) %>%
ggplot() +
geom_sf(aes(fill = factor(Cluster)), size = 0.125, color = "#b2b2b2") +
geom_sf_text(aes(label = nombre), size = 1.75, color = "black") + # Etiquetas de nombres
theme_bw() +
theme(legend.position = "right") +
scale_fill_manual(values = colores, name = "Cluster")
Modelo Normal con medidas y varianzas especificas
Una extensión natural del modelo normal multinivel consiste en permitir que tanto la media como la varianza sean específicas para cada grupo, capturando así diferencias tanto en tendencia central como en dispersión entre los grupos.
Modelo dentro de los grupos: \[ y_{i, j} \mid \theta_j, \sigma_j^2 \overset{\text{ind}}{\sim} \textsf{N}(\theta_j, \sigma_j^2) \]
Modelo entre los grupos: \[ \theta_j \mid \mu, \tau^2 \overset{\text{iid}}{\sim} \textsf{N}(\mu, \tau^2) \qquad \sigma_j^2 \mid \nu, \sigma^2 \overset{\text{iid}}{\sim} \textsf{GI}\left(\tfrac{\nu}{2}, \tfrac{\nu\,\sigma^2}{2}\right) \]
Distribución previa conjunta: \[ p(\mu, \tau^2, \nu, \sigma^2) = p(\mu)\,p(\tau^2)\,p(\nu)\,p(\sigma^2) \] con: \[ \mu \sim \textsf{N}(\mu_0, \gamma_0^2), \qquad \tau^2 \sim \textsf{GI}\left(\tfrac{\eta_0}{2}, \tfrac{\eta_0\,\tau_0^2}{2}\right), \qquad p(\nu) \propto e^{-\lambda_0 \nu}, \qquad \sigma^2 \sim \textsf{Gamma}(\alpha_0, \beta_0) \]
Los parámetros del modelo son: \[ \boldsymbol{\theta} = (\theta_1, \ldots, \theta_m, \sigma_1^2, \ldots, \sigma_m^2, \mu, \tau^2, \nu, \sigma^2) \]
Los hiperparámetros son: \[ (\mu_0, \gamma_0^2, \eta_0, \tau_0^2, \lambda_0, \alpha_0, \beta_0) \]
Entrenamiento
La estimación de los parámetros del modelo se puede hacer por medio de un muestreador de Gibbs para obtener muestras de la distribución posterior \(p(\boldsymbol{\theta}\mid\boldsymbol{y})\).
(Ejercicio.) Distribución posterior: \[ \begin{aligned} p(\boldsymbol{\theta} \mid \boldsymbol{y}) &\propto p(\boldsymbol{y} \mid \boldsymbol{\theta})\, p(\boldsymbol{\theta}) \\ &= \prod_{j=1}^m\prod_{i=1}^{n_j} \textsf{N}\left(y_{i, j} \mid \theta_{j}, \sigma_j^{2}\right) \\ &\quad\quad\times \prod_{j=1}^m \textsf{N}\left(\theta_{j} \mid \mu, \tau^{2}\right) \times \prod_{j=1}^m \textsf{GI}\left(\sigma^2_{j} \mid \tfrac{\nu}{2}, \tfrac{\nu\,\sigma^{2}}{2}\right) \\ &\quad\quad\quad\times \textsf{N}\left(\mu \mid \mu_{0}, \gamma_{0}^{2}\right) \times \textsf{GI}\left(\tau^{2} \mid \tfrac{\eta_{0}}{2}, \tfrac{\eta_{0}\,\tau_{0}^{2}}{2}\right) \\ &\quad\quad\quad\quad\times e^{-\lambda_0\,\nu} \times \textsf{G}(\sigma^2\mid\alpha_0,\beta_0) \end{aligned} \]
(Ejercicio.) Distribuciones condicionales completas: \[ \begin{aligned} \theta_{j} \mid - &\sim\textsf{N}\left(\frac{\mu / \tau^{2} + n_{j} \bar{y}_{j} / \sigma_j^{2}}{1 / \tau^{2} + n_{j} /\sigma_j^{2}}, \frac{1}{1 / \tau^{2} + n_{j} /\sigma_j^{2}}\right) \\ \sigma_{j}^{2} \,\mid\, - &\sim\textsf{GI}\left(\frac{\nu+n_{j}}{2}, \frac{\nu \sigma^{2}+\sum_{i=1}^{n_j}\left(y_{i, j}-\theta_{j}\right)^{2}}{2}\right) \\ \mu \mid - &\sim\textsf{N}\left(\frac{\mu_{0} / \gamma_{0}^{2} + m \bar{\theta} / \tau^{2}}{1 / \gamma_{0}^{2} + m / \tau^{2}}, \frac{1}{1 / \gamma_{0}^{2} + m / \tau^{2}}\right) \\ \tau^{2} \mid - &\sim\textsf{ GI }\left(\frac{\eta_{0}+m}{2}, \frac{\eta_{0} \tau_{0}^{2}+\sum_{j=1}^m\left(\theta_{j}-\mu\right)^{2}}{2}\right)\\ \sigma^{2} \mid -&\sim\textsf{G}\left(\alpha_0+\frac{m \nu}{2}, \beta_0+\frac{\nu}{2} \sum_{j=1}^m \frac{1}{\sigma_{j}^2}\right) \end{aligned} \]
La distribución condicional completa de \(\nu\) no tiene una forma cerrada conocida: \[ p\left(\nu \mid -\right) \propto\left[\frac{\left(\nu\,\sigma^{2} / 2\right)^{\nu / 2}}{\Gamma\left(\nu / 2\right)}\right]^{m}\left[\prod_{j=1}^m \frac{1}{\sigma_j^{2}}\right]^{\nu / 2} {\exp}\left\{-\nu\left(\lambda_0 + \frac{\sigma^{2}}{2} \sum_{j=1}^m \frac{1}{\sigma_{j}^{2}}\right)\right\}, \] o en escala log, \[ \log p\left(\nu \mid -\right) \propto \frac{m\,\nu}{2} \log(\nu\,\sigma^{2} / 2) - m\log\Gamma(\nu/2) -\frac{\nu}{2} \sum_{j=1}^m \log(\sigma_j^{2}) - \nu\left(\lambda_0 + \frac{\sigma^{2}}{2} \sum_{j=1}^m \frac{1}{\sigma_{j}^{2}}\right). \]
Para generar valores de \(p(\nu \mid \text{resto})\) se calcula esta distribución en escala log para un rango de valores de \(\nu\), se normaliza la distribución discreta resultante, y luego se muestrea un valor del rango de valores establecido de acuerdo con las probabilidades obtenidas.
Ejemplo: Puntajes de Matemáticas (cont.)
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), ]
# Dimensiones
dim(dat)## [1] 5511 83##
## 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))) ## [1] 33## [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)## # 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.Distribución previa
Teniendo en cuenta la información de la prueba, el modelo se ajusta utilizando los siguientes hiperparámetros: \[ \mu_0 = 50\,,\qquad \gamma^2_0 = 10^2\,,\qquad \eta_0 = 1\,,\qquad \tau^2_0 = 10^2\,,\qquad \lambda_0 = 1\,,\qquad \alpha_0 = 1\,,\qquad \beta_0 = 1/10^2\,. \]
Ajuste del modelo
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 2
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 ...Convergencia
# Cadenas: comparación de log-verosimilitud entre dos modelos
col <- RColorBrewer::brewer.pal(9, "Set1")[1:9]
# Rango común para el eje y
yrange <- range(chain1$THETA$ll, chain2$THETA$ll)
# Configuración de gráficos lado a lado
par(mfrow = c(1, 2), mar = c(2.75, 2.75, 1.5, 0.5), mgp = c(1.7, 0.7, 0))
# Modelo 1
plot(chain1$THETA$ll, type = "p", pch = ".", cex = 1.1, col = col[2],
ylim = yrange, xlab = "Iteración", ylab = "Log-verosimilitud",
main = expression(italic("Modelo 1")))
abline(h = mean(chain1$THETA$ll), lwd = 3, col = col[2])
# Modelo 2
plot(chain2$THETA$ll, type = "p", pch = ".", cex = 1.1, col = col[1],
ylim = yrange, xlab = "Iteración", ylab = "Log-verosimilitud",
main = expression(italic("Modelo 2")))
abline(h = mean(chain2$THETA$ll), lwd = 3, col = col[1])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6157 9324 10000 9521 10000 11632# Coeficiente de variación de Monte Carlo (%)
EEMC2 <- apply(X = chain2$THETA, MARGIN = 2, FUN = sd)/sqrt(neff2)
CVMC2 <- 100*abs(EEMC2/colMeans(chain2$THETA))
round(summary(CVMC2), 4)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0003 0.0212 0.0702 0.1085 0.1264 0.5342
Inferencia
Ejercicios
Considere el modelo normal jerárquico con medias específicas por grupo y varianza común. Demuestre que la varianza marginal de las observaciones está dada por \(\textsf{Var}(y_{i,j}) = \sigma^2 + \tau^2.\).
Considere el modelo normal jerárquico con medias específicas por grupo y varianza común. Utilizando los datos de la prueba Saber 11 del segundo semestre de 2023:
- Obtenga la distribución predictiva posterior para una nueva observación de en un grupo observado de su interés (e.g., Bogotá).
- Obtenga la distribución predictiva posterior para una nueva observación de en un grupo no observado (i.e., un departamento hipotético del que no se tienen observaciones).
- Evalúe la bondad de ajuste del modelo mediante un conjunto de estadísticos de prueba relevantes para un grupo observado de su interés (e.g., Bogotá), utilizando simulaciones de la distribución predictiva posterior.
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.
- 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”).
- Derive las distribuciones condicionales completas de todos los parámetros del modelo.
- Ajuste el modelo a los datos de la prueba Saber 11 del segundo semestre de 2023.
- Identifique los valores atípicos (outliers) utilizando las variables auxiliares.
- 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.
- Interprete y compare los resultados obtenidos con los del modelo normal jerárquico.
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.