Introducción

El modelo de regresión se enfoca en describir cómo el proceso generativo asociado a una variable aleatoria \(y\) varía en función de un conjunto de variables \(\boldsymbol{x} = (x_1, \ldots, x_p)\). Específicamente, un modelo de regresión define una estructura para \(p(y \mid \boldsymbol{x})\), es decir, la distribución condicional de \(y\) dado \(\boldsymbol{x}\).

La estimación de \(p(y \mid \boldsymbol{x})\) se lleva a cabo utilizando el vector de observaciones \(\boldsymbol{y} = (y_1, \ldots, y_n)\), recopilado bajo diversas condiciones representadas por \(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_n\), donde cada \(\boldsymbol{x}_i = (x_{i,1}, \ldots, x_{i,p})\) para \(i = 1, \ldots, n\).

Una posible solución a este problema es asumir que \(p(y \mid \boldsymbol{x})\) es una función suave de \(\boldsymbol{x}\), permitiendo que los valores de \(\boldsymbol{x}\) influyan en el proceso generativo de \(y\). Un modelo de regresión lineal constituye un caso particular de modelo para \(p(y \mid \boldsymbol{x})\), en el que se especifica que \(\textsf{E}(y \mid \boldsymbol{x})\) sigue una relación lineal en términos de un conjunto de parámetros \(\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)\). Esta relación se expresa como: \[ \textsf{E}(y \mid \boldsymbol{x}) = \int_{\mathcal{Y}} y\, p(y \mid \boldsymbol{x})\,\text{d}y = \sum_{k=1}^p \beta_k x_k = \boldsymbol{\beta}^{\textsf{T}} \boldsymbol{x}\,. \]

Regresión lineal Normal

El modelo de regresión lineal Normal asume que la variabilidad alrededor de \(\textsf{E}(y \mid \boldsymbol{x})\) sigue una distribución Normal, definida como: \[ y_i \mid \boldsymbol{x}_i, \boldsymbol{\beta}, \sigma^2 \stackrel{\text{iid}}{\sim} \textsf{N}(\boldsymbol{\beta}^{\textsf{T}} \boldsymbol{x}_i, \sigma^2) \qquad\Longleftrightarrow\qquad y_i = \boldsymbol{\beta}^{\textsf{T}} \boldsymbol{x}_i + \epsilon_i, \quad \epsilon_i \mid \sigma^2 \stackrel{\text{iid}}{\sim} \textsf{N}(0, \sigma^2), \] para \(i = 1, \ldots, n\).

Equivalentemente, en forma matricial: \[ \boldsymbol{y} \mid \mathbf{X}, \boldsymbol{\beta}, \sigma^2 \sim \textsf{N}_n(\mathbf{X} \boldsymbol{\beta}, \sigma^2 \mathbf{I}) \qquad\Longleftrightarrow\qquad \boldsymbol{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \mid \sigma^2 \sim \textsf{N}_n(\boldsymbol{0}, \sigma^2 \mathbf{I}), \] donde \(\mathbf{X} = [\boldsymbol{x}_1, \ldots, \boldsymbol{x}_n]^{\textsf{T}}\) es la matriz de diseño, \(\boldsymbol{\epsilon} = (\epsilon_1, \ldots, \epsilon_n)\) es el vector de términos de error, y \(\mathbf{I}\) es la matriz identidad de tamaño \(n \times n\).

(Ejercicio.) Esta formulación define completamente la distribución muestral de los datos, especificando la probabilidad conjunta de las observaciones \(\boldsymbol{y}\) dado el conjunto de predictores \(\mathbf{X}\), los parámetros \(\boldsymbol{\beta}\) y \(\sigma^2\): \[ \begin{align*} p(\boldsymbol{y} \mid \mathbf{X}, \boldsymbol{\beta}, \sigma^2) &\propto (\sigma^2)^{-n/2} \exp\left\{ -\frac{1}{2\sigma^2} (\boldsymbol{y} - \mathbf{X}\boldsymbol{\beta})^{\textsf{T}} (\boldsymbol{y} - \mathbf{X}\boldsymbol{\beta}) \right\} \\ &= (\sigma^2)^{-n/2} \exp\left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \boldsymbol{\beta}^{\textsf{T}} \boldsymbol{x}_i)^2 \right\}\,. \end{align*} \]

Mínimos cuadrados ordinarios

Dado el conjunto de datos \(\boldsymbol{y}\) y \(\mathbf{X}\), el valor óptimo de \(\boldsymbol{\beta}\) maximiza el exponente de la distribución muestral.

(Ejercicio.) Esto equivale a minimizar la suma de cuadrados residual (SSR, squares sum of residuals): \[ \textsf{SSR}(\boldsymbol{\beta}) = (\boldsymbol{y} - \mathbf{X}\boldsymbol{\beta})^{\textsf{T}} (\boldsymbol{y} - \mathbf{X}\boldsymbol{\beta}) = \boldsymbol{\beta}^{\textsf{T}} \mathbf{X}^{\textsf{T}} \mathbf{X} \boldsymbol{\beta} - 2 \boldsymbol{\beta}^{\textsf{T}} \mathbf{X}^{\textsf{T}} \boldsymbol{y} + \boldsymbol{y}^{\textsf{T}} \boldsymbol{y}. \]

(Ejercicio.) La minimización de esta expresión con respecto a \(\boldsymbol{\beta}\) conduce al valor: \[ \hat{\boldsymbol{\beta}}_{\text{ols}} = (\mathbf{X}^{\textsf{T}} \mathbf{X})^{-1} \mathbf{X}^{\textsf{T}} \boldsymbol{y}, \] conocido como el estimador de mínimos cuadrados ordinarios (OLS, ordinary least squares) de \(\boldsymbol{\beta}\).

(Ejercicio.) Adicionalmente, un estimador insesgado de \(\sigma^2\) está dado por: \[ \hat{\sigma}^2_{\text{ols}} = \frac{1}{n - p} \, \textsf{SSR}(\hat{\boldsymbol{\beta}}_{\text{ols}}), \] donde \(n\) es el número de observaciones y \(p\) es el número de parámetros estimados (incluyendo el intercepto, si corresponde).

Inferencia Bayesiana

Motivación:

  • Construcción de intervalos de credibilidad para los parámetros \(\beta_j\).
  • Cálculo de probabilidades posteriores de la forma \(\textsf{Pr}(\beta_j > 0 \mid \mathbf{X}, \boldsymbol{y})\).
  • Debido a la previa, el enfoque Bayesiano es más conservador similar a una regresión regularizada.
  • Proporciona un marco natural para la selección de modelos, al permitir comparar y evaluar modelos mediante criterios basados en la probabilidad posterior.

Modelo semi-conjugado

Distribución muestral: \[ \boldsymbol{y} \mid \mathbf{X},\boldsymbol{\beta},\sigma^2 \sim \textsf{N}_n(\mathbf{X}\boldsymbol{\beta},\sigma^2\mathbf{I}) \]

Distribución previa: \[ \boldsymbol{\beta} \sim \textsf{N}_p(\boldsymbol{\beta}_0, \mathbf{\Sigma}_0)\,,\qquad \sigma^2 \sim \textsf{GI}\left( \tfrac{\nu_0}{2}, \tfrac{\nu_0\sigma^2_0}{2} \right) \]

Los parámetros del modelo son \((\boldsymbol{\beta},\sigma^2)\).

Los hiperparámetros del modelo son \((\boldsymbol{\beta}_0, \mathbf{\Sigma}_0, \nu_0, \sigma^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{\beta},\sigma^2 \mid \mathbf{X},\boldsymbol{y})\).

  • (Ejercicio.) Distribución posterior: \[ p(\boldsymbol{\beta},\sigma^2 \mid \mathbf{X},\boldsymbol{y}) \propto \textsf{N}_n(\boldsymbol{y}\mid \mathbf{X}\boldsymbol{\beta},\sigma^2\mathbf{I}) \times \textsf{N}_p(\boldsymbol{\beta}\mid\boldsymbol{\beta}_0, \mathbf{\Sigma}_0) \times \textsf{GI}\left(\sigma^2\mid \nu_0/2, \nu_0\sigma^2_0/2 \right) \]

  • (Ejercicio.) Distribuciones condicionales completas:

\[ \begin{align*} \boldsymbol{\beta}\mid - &\sim \textsf{N}_p\left( (\mathbf{\Sigma}_0^{-1} + \tfrac{1}{\sigma^2}\mathbf{X}^{\textsf{T}}\mathbf{X} )^{-1} (\mathbf{\Sigma}_0^{-1}\boldsymbol{\beta}_0 + \tfrac{1}{\sigma^2}\mathbf{X}^{\textsf{T}}\boldsymbol{y} ) , (\mathbf{\Sigma}_0^{-1} + \tfrac{1}{\sigma^2}\mathbf{X}^{\textsf{T}}\mathbf{X} )^{-1} \right) \\ \sigma^2 \mid - &\sim \textsf{GI}\left( \frac{\nu_0 + n}{2}, \frac{\nu_0\sigma^2_0 + \textsf{SSR}(\boldsymbol{\beta})}{2} \right) \end{align*} \]

Observe que:

  • Si \(\mathbf{\Sigma}_0 \approx \mathbf{0}\,\), entonces \(\textsf{E}(\boldsymbol{\beta}\mid\mathbf{X}, \boldsymbol{y}) \approx \hat{\boldsymbol{\beta}}_{\text{ols}}\).
  • Si \(\sigma^2 \approx \infty\), entonces \(\textsf{E}(\boldsymbol{\beta}\mid\mathbf{X}, \boldsymbol{y}) \approx \boldsymbol{\beta}_0\).

Previa unitaria

Determinar valores adecuados para los hiperparámetros que representen información previa sustantiva es un desafío. Esta dificultad aumenta considerablemente a medida que crece \(p\), ya que el número de hiperparámetros se incrementa de forma drástica.

En ausencia de información previa relevante, es recomendable emplear una distribución previa que sea lo menos informativa posible.

Una opción común en estos casos es la distribución previa de información unitaria (unit information prior), que encapsula únicamente la información equivalente a la proporcionada por una sola observación.

Dado que \(\textsf{E}(\hat{\boldsymbol{\beta}}_{\text{ols}}) = \boldsymbol{\beta}\) y \((\,\textsf{Var}(\hat{\boldsymbol{\beta}}_{\text{ols}})\,)^{-1} = \frac{1}{\sigma^2} (\mathbf{X}^{\textsf{T}}\mathbf{X})\), se propone utilizar los siguientes valores: \[ \boldsymbol{\beta}_0 = \hat{\boldsymbol{\beta}}_{\text{ols}}, \qquad \mathbf{\Sigma}_0 = n\,\hat{\sigma}^2_{\text{ols}} (\mathbf{X}^{\textsf{T}} \mathbf{X})^{-1}, \qquad \nu_0 = 1, \qquad \sigma^2_0 = \hat{\sigma}^2_{\text{ols}}. \]

Estrictamente hablando, esta no es una distribución previa en el sentido convencional, ya que se basa en los datos observados. Sin embargo, dado que solo utiliza una pequeña fracción de la información disponible, equivalente a \(\frac{1}{n}\), puede interpretarse como la distribución previa de un analista con información previa débil. Este enfoque permite incorporar suavemente la información de los datos mientras se mantiene una perspectiva no informativa.

Referencias:

  • Kass, R. E., & Wasserman, L. (1995). A reference Bayesian test for nested hypotheses and its relationship to the Schwarz criterion. Journal of the american statistical association, 90(431), 928-934.
  • Kass, R. E., & Wasserman, L. (1996). The selection of prior distributions by formal rules. Journal of the American statistical Association, 91(435), 1343-1370.

Previa g

Otro principio para definir una distribución previa consiste en garantizar que la estimación de los parámetros sea invariante a cambios en la escala de los regresores. Esta propiedad se cumple siempre que \(\boldsymbol{\beta}_0 = \boldsymbol{0}\) y \(\mathbf{\Sigma}_0 = k\,(\mathbf{X}^{\textsf{T}} \mathbf{X})^{-1}\), donde \(k > 0\) actúa como un parámetro de escala que controla la dispersión de la distribución previa.

La distribución previa \(g\) (g-prior) propone que \(k = g\,\sigma^2\), con \(g > 0\). Este enfoque incorpora un factor \(g\) que controla la relación entre la varianza previa y la varianza del error. En particular, cuando \(g = n\), la distribución previa refleja la información equivalente a la proporcionada por una sola observación.

(Ejercicio.) Bajo este esquema, la distribución posterior de \(\boldsymbol{\beta}\) se expresa como: \[ \boldsymbol{\beta} \mid - \sim \textsf{N}_p\left( \frac{g}{g+1} (\mathbf{X}^{\textsf{T}} \mathbf{X})^{-1} \mathbf{X}^{\textsf{T}} \boldsymbol{y}, \frac{g}{g+1} \sigma^2 (\mathbf{X}^{\textsf{T}} \mathbf{X})^{-1} \right). \]

Este enfoque garantiza una integración coherente entre la información previa y los datos observados, ajustando la influencia de la información previa mediante el parámetro \(g\).

Bajo la previa \(g\), para simular de manera eficiente la distribución posterior de \((\boldsymbol{\beta}, \sigma^2)\), se puede aprovechar la factorización: \[ p(\boldsymbol{\beta}, \sigma^2 \mid \mathbf{X}, \boldsymbol{y}) = p(\boldsymbol{\beta} \mid \sigma^2, \mathbf{X}, \boldsymbol{y}) \, p(\sigma^2 \mid \mathbf{X}, \boldsymbol{y}). \]

La distribución marginal \(p(\sigma^2 \mid \mathbf{X}, \boldsymbol{y})\) se deriva observando que: \[ p(\sigma^2 \mid \mathbf{X}, \boldsymbol{y}) \propto p(\boldsymbol{y} \mid \sigma^2, \mathbf{X}) \, p(\sigma^2), \] donde: \[ p(\boldsymbol{y} \mid \sigma^2, \mathbf{X}) = \int_{\mathbb{R}^p} p(\boldsymbol{y} \mid \boldsymbol{\beta}, \sigma^2, \mathbf{X}) \, p(\boldsymbol{\beta} \mid \sigma^2, \mathbf{X}) \, \text{d}\boldsymbol{\beta}. \] (Ejercicio.) Al realizar esta integración, se obtiene que: \[ p(\sigma^2 \mid \mathbf{X}, \boldsymbol{y}) \propto \textsf{GI}\left(\sigma^2 \mid \frac{\nu_0 + n}{2}, \frac{\nu_0 \sigma^2_0 + \textsf{SSR}_g}{2}\right), \] donde: \[ \textsf{SSR}_g = \boldsymbol{y}^{\textsf{T}} \left(\mathbf{I}_n - \frac{g}{g+1} \mathbf{X} (\mathbf{X}^{\textsf{T}} \mathbf{X})^{-1} \mathbf{X}^{\textsf{T}}\right) \boldsymbol{y}. \]

Para simular de la distribución posterior de \((\boldsymbol{\beta}, \sigma^2)\), se sigue un sencillo algoritmo de Monte Carlo (¡no es MCMC!):

  1. Simular \(\sigma^2 \sim p(\sigma^2 \mid \mathbf{X}, \boldsymbol{y})\).
  2. Condicional en el valor simulado de \(\sigma^2\), simular \(\boldsymbol{\beta} \sim p(\boldsymbol{\beta} \mid \sigma^2, \mathbf{X}, \boldsymbol{y})\).

Este método aprovecha la conjugación inherente a la previa \(g\) para facilitar simulaciones directas de la distribución posterior.

Referencias:

  • Zellner, A. (1986). On assessing prior distributions and Bayesian regression analysis with g-prior distributions. Bayesian inference and decision techniques.

Previa difusa independiente

La previa difusa independiente es una opción estándar cuando no se dispone de información previa sustantiva o se desea minimizar su influencia en los resultados.

Este enfoque asigna distribuciones independientes a los parámetros \(\boldsymbol{\beta}\) y \(\sigma^2\) con varianzas amplias, reflejando incertidumbre sustancial. Por ejemplo, \[ \boldsymbol{\beta}_0 = \boldsymbol{0}, \quad \mathbf{\Sigma}_0 = 100\,\mathbf{I}_p, \quad \nu_0 = 1, \quad \sigma^2_0 = 100, \] siempre que 100 sea un valor “grande” respecto a la variabilidad de \(y_1,\ldots,y_n\).

Ejemplo: Absorción de oxígeno

Estudio sobre el efecto de dos regímenes de ejercicio en la absorción de oxígeno.

El estudio evaluó el efecto de dos programas de entrenamiento sobre la absorción máxima de oxígeno en un grupo de 12 hombres. Los participantes se asignaron aleatoriamente en dos grupos:

  1. Programa de carrera en terreno plano durante 12 semanas (6 participantes).
  2. Programa de aeróbicos durante 12 semanas (6 participantes).

Se midió el consumo máximo de oxígeno (en litros por minuto) de cada sujeto al correr en una cinta inclinada, tanto antes como después del programa. Este consumo se considera el mejor indicador de la capacidad de trabajo físico.

El objetivo es analizar cómo el cambio en la absorción máxima de oxígeno \(y\) depende del tipo de programa de entrenamiento.

Referencias:

  • Kuehl, R. O. (2000). Designs of experiments: Statistical principles of research design and analysis (2nd ed.). Pacific Grove, CA: Duxbury Press.
  • Hoff, P. D. (2009). A first course in Bayesian statistical methods. New York: Springer.

Variables del modelo

  • \(n\): tamaño de la muestra (12 participantes).
  • \(y\): cambio en la absorción máxima de oxígeno.
  • \(\text{trat}\): indicador del programa de entrenamiento (1 si aeróbicos, 0 si carrera).
  • \(\text{edad}\): edad del participante (en años).

El predictor lineal es de la forma: \[ \textsf{E}(y \mid \mathbf{X}) = \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \beta_4x_4, \] donde:

  • \(x_1 = 1\) (intercepto).
  • \(x_2 = 1\) si el participante está en el programa de aeróbicos, 0 si está en el programa de carrera.
  • \(x_3 = \text{edad}\).
  • \(x_4 = x_2 \cdot x_3\) (interacción entre el programa y la edad).

Coeficientes

El modelo permite estimar el cambio promedio en la absorción máxima de oxígeno según el tipo de programa y la edad del participante:

Programa de carrera: \[ \textsf{E}(y \mid \mathbf{X}) = \beta_1 + \beta_3 \cdot \text{edad}. \]

Programa de aeróbicos: \[ \textsf{E}(y \mid \mathbf{X}) = (\beta_1 + \beta_2) + (\beta_3 + \beta_4) \cdot \text{edad}. \]

# Datos
trat <- c(0, 0, 0, 0, 0, 0, 
          1, 1, 1, 1, 1, 1)
edad <- c(23, 22, 22, 25, 27, 20,
          31, 23, 27, 28, 22, 24)
y    <- c(-0.87, -10.74, -3.27, -1.97, 7.50, -7.25,
          17.05,   4.96, 10.40, 11.05, 0.26,  2.51)

# Respuesta y matriz de diseño
y <- as.matrix(y)
X <- cbind(1, trat, edad, trat*edad)
colnames(X) <- paste0("x", 1:4)

# Dimensiones
(n <- dim(X)[1])
## [1] 12
(p <- dim(X)[2])
## [1] 4

# Estimación de los coeficientes beta OLS
beta_ols <- solve(t(X) %*% X) %*% t(X) %*% y
print(round(beta_ols, 3))
##       [,1]
## x1 -51.294
## x2  13.107
## x3   2.095
## x4  -0.318
# Estimación de la varianza del error OLS
residuos   <- y - X %*% beta_ols
sig2_ols   <- sum(residuos^2) / (n - p)
print(round(sig2_ols, 3))
## [1] 8.542
# Alternativamente, se puede usar la función lm()
# fit_ols <- lm(y ~ -1 + X)  # El "-1" omite el intercepto si X ya lo incluye
# summary(fit_ols)$coefficients
# summary(fit_ols)$sigma^2
# Hiperparámetros (previa g)
nu0 <- 1
s20 <- sig2_ols
g   <- n
# Ajuste del modelo
H_g    <- (g / (g + 1)) * X %*% solve(t(X) %*% X) %*% t(X)
SSR_g  <- as.numeric(t(y) %*% (diag(n) - H_g) %*% y)
V_beta <- (g / (g + 1)) * solve(t(X) %*% X)
E_beta <- V_beta %*% t(X) %*% y

# Número de muestras
B <- 10000

# Almacenamiento
sig2_mc  <- matrix(NA, nrow = B, ncol = 1)
beta_mc  <- matrix(NA, nrow = B, ncol = p)

# Monte Carlo
set.seed(123)
for (b in 1:B) {
  sig2_mc[b]   <- 1 / rgamma(1, shape = 0.5 * (nu0 + n), rate = 0.5 * (nu0 * s20 + SSR_g))
  beta_mc[b, ] <- mvtnorm::rmvnorm(1, mean = E_beta, sigma = sig2_mc[b] * V_beta)
}

colnames(beta_mc) <- paste0("beta", 1:p)
# Estadísticos resumen para sigma^2
resumen_sigma2 <- c(
  Media = mean(sqrt(sig2_mc)),
  quantile(sqrt(sig2_mc), probs = c(0.025, 0.975))
)

# Estadísticos resumen para cada beta_j
resumen_beta <- apply(beta_mc, 2, function(x) {
  c(Media = mean(x), quantile(x, probs = c(0.025, 0.975)))
})

# Unir todo en una sola tabla
tabla_resultados <- rbind(
  sigma = round(resumen_sigma2, 3),
  round(t(resumen_beta), 3)
)

# Mostrar tabla
tabla_resultados
##         Media    2.5%   97.5%
## sigma   3.397   2.309   5.130
## beta1 -47.634 -75.525 -19.359
## beta2  12.383 -23.159  47.882
## beta3   1.946   0.725   3.142
## beta4  -0.305  -1.770   1.158
# Calculo de probabilidades sobre los coeficientes
round(colMeans(beta_mc > 0), 3)
## beta1 beta2 beta3 beta4 
## 0.001 0.768 0.999 0.330

Interpretación

Las distribuciones posteriores de \(\beta_2\) y \(\beta_4\) ofrecen evidencia limitada sobre diferencias entre los grupos, ya que sus intervalos de credibilidad incluyen el cero. La edad se asocia positivamente con el cambio en absorción en ambos programas, pero no hay evidencia sólida de que esta relación difiera según el tipo de entrenamiento.

Sin embargo, analizar estos parámetros de manera aislada no es suficiente para obtener conclusiones sólidas, ya que es necesario considerar si existen diferencias en función de la edad.

Efecto del tratamiento aeróbico \(\delta\): \[ \begin{align*} \delta(\text{edad}) &= \textsf{E}(y \mid \text{edad}, \text{aeróbico}) - \textsf{E}(y \mid \text{edad}, \text{carrera}) \\ &= ( (\beta_1 + \beta_2) + (\beta_3 + \beta_4)*\text{edad} ) - ( \beta_1 + \beta_3*\text{edad} ) \\ &= \beta_2 + \beta_4*\text{edad} \end{align*} \]

# Rango de edades observadas
rango_edad <- min(edad):max(edad)
n_edad     <- length(rango_edad)

# Inicializar matriz para almacenar muestras Monte Carlo de delta(edad)
TE <- matrix(NA, nrow = B, ncol = n_edad)

# Calcular delta(edad) = beta2 + beta4 * edad para cada edad en el rango
for (j in seq_along(rango_edad)) {
  edad_j  <- rango_edad[j]
  TE[, j] <- beta_mc[, 2] + beta_mc[, 4] * edad_j
}

Regularización

Puede haber muchas variables independientes \(x\), aunque solo algunas estén realmente asociadas con la variable dependiente \(y\). Incluir todas conduce a modelos saturados, poco parsimoniosos y de bajo rendimiento. Por ello, se recomienda conservar únicamente las variables con evidencia sólida de asociación, obteniendo así modelos más simples y con mejores propiedades de estimación y predicción.

La regresión Ridge y la regresión Lasso buscan evitar el sobreajuste, mejorar la predicción y seleccionar variables en contextos con alta dimensionalidad o multicolinealidad.

Ridge

Se basa en utilizar una distribución Normal sobre \(\boldsymbol{\beta}\), lo cual equivale a una penalización de tipo \(\ell_2\) en la perspectiva frecuentista.

Distribución muestral: \[ \boldsymbol{y} \mid \mathbf{X}, \boldsymbol{\beta}, \sigma^2 \sim \textsf{N}_n(\mathbf{X} \boldsymbol{\beta}, \sigma^2 \mathbf{I}) \]

Previa sobre los coeficientes de regresión: \[ \beta_j \mid \lambda^2 \sim \textsf{N}\left(0, \frac{1}{\lambda^2}\right), \qquad j = 1, \dots, p \]

Previa sobre la varianza del error: \[ \sigma^2 \sim \textsf{IG}(a_\sigma, b_\sigma) \]

Previa sobre el parámetro de regularización: \[ \lambda^2 \sim \textsf{Gamma}(a_\lambda, b_\lambda) \]

El parámetro \(\lambda\) controla el grado de regularización: valores grandes implican mayor penalización (más shrinkage hacia cero).

(Ejercicio.) El máximo a posteriori (MAP) coincide con la solución de Ridge en su forma frecuentista. La distribución condicional completa (en escala log) de \(\boldsymbol{\beta}\) es proporcial a: \[ \log p(\boldsymbol{\beta} \mid -) \propto -\frac{1}{2\sigma^2} \| \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \|^2 - \frac{\lambda^2}{2} \| \boldsymbol{\beta} \|^2, \] por lo que maximizar esta expresión (calcular el MAP) equivale a minimizar \[ \| \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \|^2 + \lambda_{\text{ridge}} \| \boldsymbol{\beta} \|^2 \] donde \(\lambda_{\text{ridge}} = \sigma^2\lambda^2\). Esta es exactamente la forma de la regresión Ridge frecuentista, por lo que el MAP coincide con la solución de Ridge.

Distribuciones condicionales completas

(Ejercicio.) Distribuciones condicionales completas:

  • \(\boldsymbol{\beta} \mid - \sim \textsf{N}_p(\mathbf{m}, \mathbf{V})\), donde: \[ \mathbf{m} = \left( \mathbf{X}^\top \mathbf{X} + \sigma^2 \lambda^2 \mathbf{I} \right)^{-1} \mathbf{X}^\top \boldsymbol{y} , \quad \mathbf{V} = \sigma^2\left( \mathbf{X}^\top \mathbf{X} + \sigma^2 \lambda^2 \mathbf{I} \right)^{-1}\,. \]

  • \(\sigma^2 \mid - \sim \textsf{IG}(\text{A}, \text{B})\), donde: \[ \text{A} = a_\sigma + \frac{n}{2}, \quad \text{B} = b_\sigma + \frac{ \| \boldsymbol{y} - \mathbf{X} \boldsymbol{\beta} \|^2 }{2}\,. \]

  • \(\lambda^2 \mid - \sim \textsf{Gamma}(\text{A}, \text{B})\), donde: \[ \text{A} = a_\lambda + \frac{p}{2},\qquad \text{B} = b_\lambda + \frac{\| \boldsymbol{\beta} \|^2}{2}\,. \]

Lasso

Se basa en utilizar una distribución Laplace sobre cada coeficiente \(\beta_j\), lo cual equivale a una penalización de tipo \(\ell_1\) en la perspectiva frecuentista:

Distribución muestral: \[ \boldsymbol{y} \mid \mathbf{X}, \boldsymbol{\beta}, \sigma^2 \sim \textsf{N}_n(\mathbf{X} \boldsymbol{\beta}, \sigma^2 \mathbf{I}) \]

Previa sobre los coeficientes de regresión: \[ \beta_j\mid\lambda \sim \textsf{Laplace}\left(0,\frac{1}{\lambda}\right), \qquad j = 1, \dots, p \] con \(p(\beta_j \mid \lambda) \propto \exp\left(-\lambda|\beta_j|\right)\).

(Ejercicio.) Esto es equivalente a \[ \beta_j \mid \tau_j^2 \sim \textsf{N}(0, \tau_j^2), \qquad \tau_j^2\mid \lambda^2 \sim \textsf{Exp}\left( \frac{\lambda^2}{2} \right), \qquad j = 1, \dots, p \]

Previa sobre la varianza del error: \[ \sigma^2 \sim \textsf{GI}(a_\sigma, b_\sigma) \]

Previa sobre el parámetro de regularización: \[ \lambda^2 \sim \textsf{Gamma}(a_\lambda, b_\lambda) \]

La penalización del Lasso es más agresiva que la del Ridge porque la distribución Laplace tiene colas más pesadas y un pico más agudo en cero que la distribución Normal.

(Ejercicio.) El MAP coincide con la solución de Lasso en su forma frecuentista. La distribución condicional completa (en escala log) de \(\boldsymbol{\beta}\) es proporcional a: \[ \log p(\boldsymbol{\beta} \mid -) \propto -\frac{1}{2\sigma^2} \| \boldsymbol{y} - \mathbf{X} \boldsymbol{\beta} \|^2 - \lambda \sum_{j=1}^p |\beta_j|, \] por lo que maximizar esta expresión equivale a minimizar \[ \| \boldsymbol{y} - \mathbf{X} \boldsymbol{\beta} \|^2 + \lambda_{\text{lasso}} \| \boldsymbol{\beta} \|_1, \] donde \(\lambda_{\text{lasso}} = 2\sigma^2\lambda\). Esta es exactamente la forma de la regresión Lasso frecuentista, por lo que el MAP coincide con la solución de Lasso.

Distribuciones condicionales completas

(Ejercicio.) Distribuciones condicionales completas:

  • \(\boldsymbol{\beta} \mid - \sim \textsf{N}_p(\mathbf{m}, \mathbf{V})\), donde: \[ \mathbf{m} = \left( \mathbf{X}^\top \mathbf{X} + \sigma^2\mathbf{D}_\tau^{-1} \right)^{-1} \mathbf{X}^\top \boldsymbol{y}, \qquad \mathbf{V} = \sigma^2\left( \mathbf{X}^\top \mathbf{X} + \sigma^2\mathbf{D}_\tau^{-1} \right)^{-1}\,. \] y \(\mathbf{D}_{\tau} = \text{diag}(\tau_1^2, \dots, \tau_p^2)\).

  • \(\tau_j^2 \mid - \sim \textsf{Gaussiana-Inversa-Generalizada}(\nu, \text{A}, \text{B})\), donde: \[ \nu = \frac12,\qquad \text{A} = \lambda^2,\qquad \text{B} = \beta^2_j\,. \]

  • \(\sigma^2 \mid - \textsf{GI}(\text{A}, \text{B})\), donde: \[ \text{A} = a_\sigma + \frac{n}{2},\qquad \text{B} = \ b_\sigma + \frac{\| \boldsymbol{y} - \mathbf{X} \boldsymbol{\beta} \|^2}{2} \,. \]

  • \(\lambda^2 \mid - \textsf{Gamma}(\text{A}, \text{B})\), donde: \[ \text{A} = a_\lambda + p,\qquad \text{B} =\ b_\lambda + \frac{1}{2} \sum_{j=1}^p \tau_j^2\,. \]

Ejemplo: Simulación de regularización

Los datos se generan a partir del modelo de regresión lineal: \[ \boldsymbol{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim \textsf{N}_n(\boldsymbol{0}, \sigma^2 \mathbf{I}) \] donde \(\mathbf{X} \in \mathbb{R}^{n \times p}\) es la matriz de predictores, \(\boldsymbol{\beta} \in \mathbb{R}^p\) es el vector de coeficientes verdaderos, y \(\boldsymbol{y} \in \mathbb{R}^n\) es la variable respuesta.

Se simulan \(n = 50\) observaciones y \(p = 40\) predictores independientes estándar. El vector \(\boldsymbol{\beta}\) contiene 5 coeficientes no nulos \((3, -2, 1.5, -1, 2)\) y 15 coeficientes exactamente iguales a cero. El término de error sigue una distribución normal con varianza \(\sigma^2 = 3\).

Este diseño permite evaluar el desempeño de los modelos Bayesianos Ridge y Lasso en presencia de sparsity.

Simulación de datos

# Semilla para reproducibilidad
set.seed(123)

# Parámetros del diseño
n <- 50 
p <- 40 
s <- 5 

# Matriz de diseño
X <- matrix(rnorm(n * p), nrow = n, ncol = p)

# Coeficientes verdaderos: sparse
beta_true <- c(3, -2, 1.5, -1, 2, rep(0, p - s))

# Ruido
sigma2_true <- 3
eps_true    <- rnorm(n, mean = 0, sd = sqrt(sigma2_true))

# Variable respuesta
y <- X %*% beta_true + eps_true

Entrenamiento: Previa g

gibbs_gprior <- function(y, X, B, g, nu0, s20) {
  # Dimensiones
  n <- length(y)
  p <- ncol(X)

  # Constantes
  H_g    <- (g / (g + 1)) * X %*% solve(t(X) %*% X) %*% t(X)
  SSR_g  <- as.numeric(t(y) %*% (diag(n) - H_g) %*% y)
  V_beta <- (g / (g + 1)) * solve(t(X) %*% X)
  E_beta <- V_beta %*% t(X) %*% y

  # Almacenamiento
  beta_samples   <- matrix(NA, nrow = B, ncol = p)
  sigma2_samples <- numeric(B)
  loglik_samples <- numeric(B)

  # Muestreo directo
  for (b in 1:B) {
    sigma2 <- 1 / rgamma(1, shape = 0.5 * (nu0 + n), rate = 0.5 * (nu0 * s20 + SSR_g))
    beta   <- c(mvtnorm::rmvnorm(1, mean = E_beta, sigma = sigma2 * V_beta))
    
    # Almacenamiento
    beta_samples  [b, ] <- beta
    sigma2_samples[b]   <- sigma2
    loglik_samples[b]   <- sum(dnorm(y, mean = c(X %*% beta), sd = sqrt(sigma2), log = TRUE))
    
    # Avance cada 10%
    if (b %% ceiling(B / 10) == 0)
      cat(paste0("Iteración ", b, " de ", B, " (", round(100 * b / B), "%)\n"))
  }

  # Etiquetas
  colnames(beta_samples) <- paste0("beta", 1:p)

  # Resultado
  return(list(
    beta   = beta_samples,
    sigma2 = sigma2_samples,
    loglik = loglik_samples
  ))
}

Entrenamiento: Ridge

gibbs_ridge <- function(y, X, B, n_burn, a_sigma, b_sigma, a_lambda, b_lambda) {
  # Dimensiones
  n <- length(y)
  p <- ncol(X)
  total_iter <- B + n_burn
  
  # Constantes
  XtX <- t(X) %*% X
  Xty <- t(X) %*% y
  
  # Almacenamiento posterior
  beta_samples    <- matrix(NA, nrow = B, ncol = p)
  sigma2_samples  <- numeric(B)
  lambda2_samples <- numeric(B)
  loglik_samples  <- numeric(B)
  
  # Inicialización
  beta    <- rep(0, p)
  sigma2  <- 1
  lambda2 <- 1
  
  # Iteraciones del Gibbs sampler
  for (b in 1:total_iter) {
    # actualizar beta
    VV   <- solve(XtX + diag(sigma2 * lambda2, p))
    mm   <- VV %*% Xty
    beta <- c(mvtnorm::rmvnorm(1, mean = mm, sigma = sigma2 * VV))
       
    # actualizar sigma2
    res    <- y - X %*% beta
    AA     <- a_sigma + 0.5 * n
    BB     <- b_sigma + 0.5 * sum(res^2)
    sigma2 <- 1 / rgamma(1, shape = AA, rate = BB)
    
    # actualizar lambda2
    AA     <- a_lambda + 0.5 * p
    BB     <- b_lambda + 0.5 * sum(beta^2)
    lambda <- rgamma(1, shape = AA, rate = BB)
    
    # Guardar muestras (después del burn-in)
    if (b > n_burn) {
      i <- b - n_burn
      beta_samples   [i, ] <- beta
      sigma2_samples [i]   <- sigma2
      lambda2_samples[i]   <- lambda2
      
      # Log-verosimilitud
      loglik_samples[i] <- sum(dnorm(y, mean = c(X %*% beta), sd = sqrt(sigma2), log = TRUE))
    }
    
    # Avance cada 10%
    if (b %% ceiling(total_iter / 10) == 0)
      cat(paste0("Iteración ", b, " de ", total_iter, " (", round(100 * b / total_iter), "%)\n"))
  }
  
  # Etiquetas
  colnames(beta_samples) <- paste0("beta", 1:p)
  
  # Resultado
  return(list(
    beta    = beta_samples,
    sigma2  = sigma2_samples,
    lambda2 = lambda2_samples,
    loglik  = loglik_samples
  ))
}

Entrenamiento: Lasso

gibbs_lasso <- function(y, X, B, n_burn, a_sigma, b_sigma, a_lambda, b_lambda) {
  # Dimensiones
  n <- length(y)
  p <- ncol(X)
  total_iter <- B + n_burn
  
  # Constantes
  XtX <- t(X) %*% X
  Xty <- t(X) %*% y
  
  # Almacenamiento posterior
  beta_samples    <- matrix(NA, nrow = B, ncol = p)
  sigma2_samples  <- numeric(B)
  lambda2_samples <- numeric(B)
  loglik_samples  <- numeric(B)
  
  # Inicialización
  beta    <- rep(0, p)
  sigma2  <- 1
  lambda2 <- 1
  tau2    <- rep(1, p)
  
  # Gibbs sampler
  for (b in 1:total_iter) {
    # Actualizar beta
    VV   <- solve(XtX + diag(sigma2 / tau2))
    mm   <- VV %*% Xty
    beta <- c(mvtnorm::rmvnorm(1, mean = mm, sigma = sigma2 * VV))
    
    # Actualizar tau2_j | beta_j, lambda2
    tau2 <- GIGrvg::rgig(n = p, lambda = 0.5, chi = beta^2, psi = lambda2)
    
    # Actualizar sigma2
    AA     <- a_sigma + 0.5 * n
    BB     <- b_sigma + 0.5 * sum((y - X %*% beta)^2)
    sigma2 <- 1 / rgamma(1, shape = AA, rate = BB)
    
    # 4. Actualizar lambda2
    AA      <- a_lambda + p
    BB      <- b_lambda + 0.5 * sum(tau2)
    lambda2 <- rgamma(1, shape = AA, rate = BB)
    
    # Guardar muestras después del burn-in
    if (b > n_burn) {
      i <- b - n_burn
      beta_samples   [i, ] <- beta
      sigma2_samples [i]   <- sigma2
      lambda2_samples[i]   <- lambda2
      
      # Log-verosimilitud
      loglik_samples[i] <- sum(dnorm(y, mean = c(X %*% beta), sd = sqrt(sigma2), log = TRUE))
    }
    
    # Avance
    if (b %% ceiling(total_iter / 10) == 0) {
      cat(paste0("Iteración ", b, " de ", total_iter, " (", round(100 * b / total_iter), "%)\n"))
    }
  }
  
  colnames(beta_samples) <- paste0("beta", 1:p)
  
  return(list(
    beta    = beta_samples,
    sigma2  = sigma2_samples,
    lambda2 = lambda2_samples,
    loglik  = loglik_samples
  ))
}

Asjute de los modelos

# Estimación de los coeficientes beta OLS
beta_ols <- solve(t(X) %*% X) %*% t(X) %*% y

# Estimación de la varianza del error OLS
residuos   <- y - X %*% beta_ols
sig2_ols   <- sum(residuos^2) / (n - p)
# Reproducibilidad
set.seed(123)

# Hiperparámetros previa g
nu0 <- 1
s20 <- sig2_ols
g   <- n

# Ajuste del modelo
B <- 10000

chain_gprior <- gibbs_gprior(y, X, B, g, nu0, s20) 
## Iteración 1000 de 10000 (10%)
## Iteración 2000 de 10000 (20%)
## Iteración 3000 de 10000 (30%)
## Iteración 4000 de 10000 (40%)
## Iteración 5000 de 10000 (50%)
## Iteración 6000 de 10000 (60%)
## Iteración 7000 de 10000 (70%)
## Iteración 8000 de 10000 (80%)
## Iteración 9000 de 10000 (90%)
## Iteración 10000 de 10000 (100%)
# Reproducibilidad
set.seed(123)

# Hiperparámetros Ridge
a_sigma  <- 3
b_sigma  <- (a_sigma) * sig2_ols
a_lambda <- 0.001
b_lambda <- 0.001

# Ajuste del modelo
B      <- 10000
n_burn <- 1000

chain_ridge <- gibbs_ridge(y, X, B, n_burn, a_sigma, b_sigma, a_lambda, b_lambda)
## Iteración 1100 de 11000 (10%)
## Iteración 2200 de 11000 (20%)
## Iteración 3300 de 11000 (30%)
## Iteración 4400 de 11000 (40%)
## Iteración 5500 de 11000 (50%)
## Iteración 6600 de 11000 (60%)
## Iteración 7700 de 11000 (70%)
## Iteración 8800 de 11000 (80%)
## Iteración 9900 de 11000 (90%)
## Iteración 11000 de 11000 (100%)
# Reproducibilidad
set.seed(123)

# Hiperparámetros Lasso
a_sigma  <- 3
b_sigma  <- (a_sigma) * sig2_ols
a_lambda <- 0.001
b_lambda <- 0.001

# Ajuste del modelo
B      <- 10000
n_burn <- 1000

chain_lasso <- gibbs_lasso(y, X, B, n_burn, a_sigma, b_sigma, a_lambda, b_lambda)
## Iteración 1100 de 11000 (10%)
## Iteración 2200 de 11000 (20%)
## Iteración 3300 de 11000 (30%)
## Iteración 4400 de 11000 (40%)
## Iteración 5500 de 11000 (50%)
## Iteración 6600 de 11000 (60%)
## Iteración 7700 de 11000 (70%)
## Iteración 8800 de 11000 (80%)
## Iteración 9900 de 11000 (90%)
## Iteración 11000 de 11000 (100%)

Convergencia

# Log-verosimilitud real
ll_true <- sum(dnorm(y, mean = c(X %*% beta_true), sd = sqrt(sigma2_true), log = TRUE))

# Recorrido de valores de la log-verosimilitud
y_range <- range(ll_true, chain_gprior$loglik, chain_ridge$loglik, chain_ridge$loglik)

# Configuración de gráficos lado a lado
par(mfrow = c(3, 1), mar = c(2.75, 2.75, 1.5, 0.5), mgp = c(1.7, 0.7, 0))

# g prior
plot(chain_gprior$loglik, type = "p", pch = ".", cex = 1.1,
     ylim = y_range, xlab = "Iteración", ylab = "Log-verosimilitud",
     main = expression(italic("g prior")))
abline(h = ll_true, lwd = 3, col = "red")

# Ridge
plot(chain_ridge$loglik, type = "p", pch = ".", cex = 1.1,
     ylim = y_range, xlab = "Iteración", ylab = "Log-verosimilitud",
     main = expression(italic("Ridge")))
abline(h = ll_true, lwd = 3, col = "red")

# Lasso
plot(chain_lasso$loglik, type = "p", pch = ".", cex = 1.1,
     ylim = y_range, xlab = "Iteración", ylab = "Log-verosimilitud",
     main = expression(italic("Lasso")))
abline(h = ll_true, lwd = 3, col = "red")

Inferencia

# Cálculos individuales
res_gprior <- c(mean(chain_gprior$sigma2), quantile(chain_gprior$sigma2, probs = c(0.025, 0.975)))
res_ridge  <- c(mean(chain_ridge$sigma2),  quantile(chain_ridge$sigma2,  probs = c(0.025, 0.975)))
res_lasso  <- c(mean(chain_lasso$sigma2),  quantile(chain_lasso$sigma2,  probs = c(0.025, 0.975)))

# Unir en matriz
tabla_sigma2 <- rbind(
  gprior = res_gprior,
  Ridge  = res_ridge,
  Lasso  = res_lasso
)

# Asignar nombres de columnas
colnames(tabla_sigma2) <- c("Media", "IC 2.5%", "IC 97.5%")

# Mostrar tabla formateada
knitr::kable(tabla_sigma2, digits = 3, align = "c", caption = "Resumen posterior de $\\sigma^2$")
Resumen posterior de \(\sigma^2\)
Media IC 2.5% IC 97.5%
gprior 0.736 0.497 1.076
Ridge 1.709 0.869 3.311
Lasso 1.752 0.890 3.412
# Muestras
beta_gprior <- chain_gprior$beta
beta_ridge  <- chain_ridge$beta
beta_lasso  <- chain_lasso$beta

# Dimensión
p <- ncol(beta_gprior)
x_pos <- 1:p

# Estadísticas posteriores
media_gprior <- colMeans(beta_gprior)
ic95_gprior  <- apply(beta_gprior, 2, quantile, probs = c(0.025, 0.975))

media_ridge  <- colMeans(beta_ridge)
ic95_ridge   <- apply(beta_ridge, 2, quantile, probs = c(0.025, 0.975))

media_lasso  <- colMeans(beta_lasso)
ic95_lasso   <- apply(beta_lasso, 2, quantile, probs = c(0.025, 0.975))

# Rango común para los ejes y
ylim_all <- range(c(ic95_gprior, ic95_ridge, ic95_lasso, beta_true))

# Gráfico lado a lado
par(mfrow = c(3, 1), mar = c(2.75, 2.75, 1.5, 0.5), mgp = c(1.7, 0.7, 0))

# g prior
plot(NA, NA, xlab = "Coeficientes", ylab = expression(beta[j]), 
     xlim = c(0.5, p + 0.5), ylim = ylim_all, main = "g prior")

abline(h = 0, col = "gray", lwd = 2)

for (j in x_pos) {
  # IC 95%
  segments(x0 = j, y0 = ic95_gprior[1, j], x1 = j, y1 = ic95_gprior[2, j], lwd = 1)
  # Media posterior
  points(x = j, y = media_gprior[j], pch = 16, cex = 0.6)
  # Valor verdadero
  points(x = j, y = beta_true[j], pch = 15, col = "red", cex = 0.8)
}

# RIDGE
plot(NA, NA, xlab = "Coeficientes", ylab = expression(beta[j]), 
     xlim = c(0.5, p + 0.5), ylim = ylim_all, main = "Ridge")

abline(h = 0, col = "gray", lwd = 2)

for (j in x_pos) {
  # IC 95%
  segments(x0 = j, y0 = ic95_ridge[1, j], x1 = j, y1 = ic95_ridge[2, j], lwd = 1)
  # Media posterior
  points(x = j, y = media_ridge[j], pch = 16, cex = 0.6)
  # Valor verdadero
  points(x = j, y = beta_true[j], pch = 15, col = "red", cex = 0.8)
}

# LASSO
plot(NA, NA, xlab = "Coeficientes", ylab = expression(beta[j]), 
     xlim = c(0.5, p + 0.5), ylim = ylim_all, main = "Lasso")

abline(h = 0, col = "gray", lwd = 2)

for (j in x_pos) {
  # IC 95%
  segments(x0 = j, y0 = ic95_lasso[1, j], x1 = j, y1 = ic95_lasso[2, j], lwd = 1)
  # Media posterior
  points(x = j, y = media_lasso[j], pch = 16, cex = 0.6)
  # Valor verdadero
  points(x = j, y = beta_true[j], pch = 15, col = "red", cex = 0.8)
}

Bondad de ajuste

Dados los valores reales \(\beta_j\) y los valores estimados \(\hat{\beta}_j\) de los coeficientes de la regresión, para \(j = 1, \dots, p\), las métricas más comunes para comparar los modelos en términos de estimabilidad son:

  • Raíz del error cuadrático medio: \[ \text{RMSE} = \sqrt{\frac{1}{p} \sum_{j=1}^p (\beta_j - \hat{\beta}_j)^2}\,. \]

  • Error absoluto medio: \[ \text{MAE} = \frac{1}{p} \sum_{j=1}^p \left| \beta_j - \hat{\beta}_j \right|\,. \]

# OLS
beta_hat_ols <- solve(t(X) %*% X) %*% t(X) %*% y
RMSE_ols     <- sqrt(mean((beta_true - beta_hat_ols)^2))
MAE_ols      <- mean(abs(beta_true - beta_hat_ols))

# g-prior
beta_hat_gprior <- colMeans(chain_gprior$beta)
RMSE_gprior     <- sqrt(mean((beta_true - beta_hat_gprior)^2))
MAE_gprior      <- mean(abs(beta_true - beta_hat_gprior))

# Ridge
beta_hat_ridge <- colMeans(chain_ridge$beta)
RMSE_ridge     <- sqrt(mean((beta_true - beta_hat_ridge)^2))
MAE_ridge      <- mean(abs(beta_true - beta_hat_ridge))

# Lasso
beta_hat_lasso <- colMeans(chain_lasso$beta)
RMSE_lasso     <- sqrt(mean((beta_true - beta_hat_lasso)^2))
MAE_lasso      <- mean(abs(beta_true - beta_hat_lasso))

# Tabla
tab <- rbind(
  c(RMSE_ols,    MAE_ols),
  c(RMSE_gprior, MAE_gprior),
  c(RMSE_ridge,  MAE_ridge),
  c(RMSE_lasso,  MAE_lasso)
)

colnames(tab) <- c("RMSE", "MAE")
rownames(tab) <- c("OLS", "g-prior", "Ridge", "Lasso")

knitr::kable(tab, digits = 3, align = "c", caption = "Comparación de los modelos en términos de estimabilidad.")
Comparación de los modelos en términos de estimabilidad.
RMSE MAE
OLS 0.458 0.377
g-prior 0.449 0.371
Ridge 0.341 0.289
Lasso 0.427 0.353

Ejemplo: Temperatura de la superficie del mar

Se considera la media de las observaciones registradas por medio de varios dispositivos acerca de la temperatura de la superficie del mar (SST, surface sea temperature) en \(^\text{o}\)C, junto con el tipo de dispositivo (cuatro categorías: bucket, eri, d.buoy, f.buoy) con sus respectivas ubicaciones (latitud y longitud) en el mar Mediterráneo en diciembre de 2003:

  • bucket (baldes): Mediciones con baldes (buckets) lanzados desde barcos, donde la temperatura se mide con un termómetro tras recoger agua superficial. Pueden presentar sesgos por enfriamiento durante la recuperación.

  • eri (engine room intake): Mediciones con sensores en la entrada de agua a la sala de máquinas, a pocos metros de profundidad. Pueden estar afectadas por el calor del motor o la ubicación de la toma.

  • d.buoy (drifting buoy): Mediciones con boyas derivantes que flotan con las corrientes, equipadas con sensores que miden la SST y transmiten los datos vía satélite. Ofrecen amplia cobertura espacial y son clave para el monitoreo global.

  • f.buoy (fixed buoy): Mediciones con boyas fijas ancladas en un punto, que registran la SST y otras variables oceánicas y meteorológicas. Permiten observaciones continuas en ubicaciones estratégicas.

Estas categorías representan fuentes heterogéneas de medición que pueden diferir en precisión, profundidad de muestreo y sesgos sistemáticos, lo cual es relevante para análisis estadísticos que buscan modelar o corregir la SST observada.

Tratamiento de datos

La base de datos sst.dat contiene \(m = 86\) registros de la forma \[ \left[\,\text{id}_j,\text{latitud}_j,\text{longitud}_j,\bar{y}_j,n_j,\text{tipo}_j\,\right] \] donde \(\bar{y}_j=\frac{1}{n_j}\sum_{i=1}^{n_j} y_{i,j}\) es el promedio de las \(n_j\) temperaturas registradas por el dispositivo \(j\), para \(j=1,\ldots,m\).

# Datos
df <- read.table(file = "sst.dat", header = T)

# Dimensiones
dim(df)
## [1] 86  6
# Encabezado de la base de datos
head(df, n = 10)
##    id  lat  lon temp  n   type
## 1   1 33.8 33.1 21.3  1 bucket
## 2   2 34.5 25.8 19.1  2 d.buoy
## 3   3 33.9 29.3 19.9  1 bucket
## 4   4 33.3 32.9 20.6  2 d.buoy
## 5   5 34.3 33.4 20.0 23 f.buoy
## 6   6 34.1 32.2 20.5 62 f.buoy
## 7   7 32.2 32.7 21.5 37 f.buoy
## 8   8 32.4 32.8 19.2 94 f.buoy
## 9   9 32.7 32.2 21.0 24 f.buoy
## 10 10 33.5 25.3 17.6  1 bucket
# Separar datos por tipo de dispositivo
df1 <- df[df$type == "bucket", ]  # 1 = bucket
df2 <- df[df$type == "eri",    ]  # 2 = eri
df3 <- df[df$type == "d.buoy", ]  # 3 = d.buoy
df4 <- df[df$type == "f.buoy", ]  # 4 = f.buoy

# Combinar los subconjuntos en un único data frame (ordenado por dispositivo)
df <- rbind(df1, df2, df3, df4)

# Tamaños de cada subconjunto
m1 <- nrow(df1)
m2 <- nrow(df2)
m3 <- nrow(df3)
m4 <- nrow(df4)

# Codificar tipo de dispositivo como variable numérica
df$cod <- c(rep(1, m1), rep(2, m2), rep(3, m3), rep(4, m4))

# Tabla de frecuencias por tipo de dispositivo
table(factor(df$type, levels = c("bucket", "eri", "d.buoy", "f.buoy")))
## 
## bucket    eri d.buoy f.buoy 
##     36     35     10      5
# Mapa usando RgoogleMaps
suppressMessages(suppressWarnings(library(RgoogleMaps)))

center <- c(median(df$lat), median(df$lon))
zoom   <- min(MaxZoom(range(df$lat), range(df$lon)))

mymap  <- GetMap(
  center   = center,
  zoom     = zoom,
  destfile = "medsea.png",
  maptype  = "mobile"
)

PlotOnStaticMap(
  mymap,
  lat = df$lat,
  lon = df$lon,
  col = df$cod,
  pch = 20
)

Los tipos de dispositivos son \(\color{black}{\text{bucket}}\), \(\color{red}{\text{eri}}\), \(\color{green}{\text{d.buoy}}\), \(\color{blue}{\text{d.buoy}}\).

Análisis exploratorio de datos

# Número de dispositivos (grupos)
(m <- nrow(df))
## [1] 86
# Variable respuesta (SST)
y <- df$temp
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.00   18.70   19.65   19.56   20.60   22.00
# Tamaños de muestra
nj <- df$n
table(nj)
## nj
##  1  2  3 23 24 37 62 94 
## 70  7  4  1  1  1  1  1

Modelo

Los \(m\) registros se asumen como grupos independientes, cada uno con \(n_j\) observaciones \(y_{1,j},\ldots,y_{n_j,j}\) condicionalmente independientes y normalmente distribuidas con media \(\mu_j\) y varianza \(\sigma^2_j\), i.e., \(y_{i,j}\mid\mu_j,\sigma^2_j \stackrel{\text{iid}}{\sim} \textsf{N}(\mu_j,\sigma^2_j)\), lo que significa que \[ \bar{y}_j\mid\mu_j,\sigma^2_j \stackrel{\text{ind}}{\sim} \textsf{N}(\mu_j,\sigma^2_j/n_j)\,, \] donde \(\bar{y}_j=\frac{1}{n_j}\sum_{i=1}^{n_j} y_{i,j}\), para \(j = 1,\ldots,m\).

Para introducir la asociación entre la repuesta media de las observaciones \(\mu_j\) y las covariables \(x_{j,1},\ldots,x_{j,p}\), se considera el predictor lineal \(\eta_j=\boldsymbol{x}_j^{\textsf{T}}\boldsymbol{\beta} = \sum_{k=1}^p \beta_k x_{j,k}\) por medio de una dependencia estocástica Normal de la forma \[ \mu_j\mid\boldsymbol{\beta},\tau^2 \stackrel{\text{ind}}{\sim} \textsf{N}(\boldsymbol{x}_j^{\textsf{T}}\boldsymbol{\beta},\tau^2) \] donde \(\boldsymbol{x}_j=(x_{j,1},\ldots,x_{j,p})\) y \(\boldsymbol{\beta}=(\beta_1,\ldots,\beta_p)\). Esta formulación permite regularización o suavizamiento entre grupos, especialmente útil cuando algunos tienen pocas observaciones.

(Ejercicio.) Bajo esta formulación se tiene que \[ p(\bar{y}_j\mid\boldsymbol{\beta},\sigma^2_j,\tau^2) = \int_{\mathbb{R}} p(\bar{y}_j,\mu_j\mid\boldsymbol{\beta},\sigma^2_j,\tau^2)\,\text{d}\mu_j = \textsf{N}(\bar{y}_j\mid\boldsymbol{x}_j^{\textsf{T}}\boldsymbol{\beta},\sigma^2_j/n_j+\tau^2)\,. \]

En resumen, \[ \bar{\boldsymbol{y}}\mid\boldsymbol{\mu},\boldsymbol{\sigma}^2\sim\textsf{N}_m(\boldsymbol{\mu},\textsf{diag}(\sigma^2_1/n_1,\ldots,\sigma^2_m/n_m)) \qquad\text{y}\qquad \boldsymbol{\mu}\mid\boldsymbol{\beta},\tau^2\sim\textsf{N}_m(\mathbf{X}\boldsymbol{\beta},\tau^2\mathbf{I})\,, \] donde \(\bar{\boldsymbol{y}}=(\bar{y}_1,\ldots,\bar{y}_m)\), \(\boldsymbol{\mu}=(\mu_1,\ldots,\mu_m)\), \(\boldsymbol{\sigma}^2=(\sigma^2_1,\ldots,\sigma^2_m)\) y \(\mathbf{X}=(\boldsymbol{x}_1^{\textsf{T}},\ldots,\boldsymbol{x}_m^{\textsf{T}})\).

(Ejercicio.) Para la estructura de regresión se considera la previa impropia \[ p(\boldsymbol{\beta},\log\tau^2) \propto 1 \qquad\Longleftrightarrow \qquad p(\boldsymbol{\beta},\tau^2) \propto \frac{1}{\tau^2} \] la cual no requiere la especificación de hiperparámetros.

Los componentes de intra-varianza se modelan jerárquicamente mediante \[ \sigma^2_j\mid\nu,\upsilon^2 \stackrel{\text{iid}}{\sim} \textsf{GI}(\nu/2,\nu\upsilon^2/2)\,,\qquad p(\nu)\propto e^{-\lambda_0\nu} \qquad\text{y}\qquad \upsilon^2\sim\textsf{G}(\alpha_0/2,\beta_0/2)\,. \]

Los parámetros del modelo son \(\boldsymbol{\mu}\), \(\boldsymbol{\sigma}^2\), \(\boldsymbol{\beta}\), \(\tau^2\), \(\nu\) y \(\upsilon^2\).

Los hiperparámetros del modelo son \(\lambda_0\), \(\alpha_0\) y \(\beta_0\).

Distribución posterior

(Ejercicio.) La distribución posterior de \(\mathbf{\Theta} = (\boldsymbol{\mu},\boldsymbol{\sigma}^2,\boldsymbol{\beta},\tau^2,\nu,\upsilon^2)\) es: \[ \begin{align*} p(\mathbf{\Theta}\mid \bar{\boldsymbol{y}}) &\propto \prod_{j=1}^m \textsf{N}(\bar{y}_j\mid\mu_j,\sigma^2_j/n_j) \\ &\hspace{1cm}\times \prod_{j=1}^m \textsf{N}(\mu_j\mid\boldsymbol{x}_j^{\textsf{T}}\boldsymbol{\beta},\tau^2) \times \prod_{j=1}^m \textsf{GI}(\sigma^2_j\mid\nu/2,\nu\upsilon^2/2)\\ &\hspace{2cm}\times \frac{1}{\tau^2} \times e^{-\lambda_0\nu} \times \textsf{G}(\upsilon^2\mid\alpha_0/2,\beta_0/2)\,. \end{align*} \]

Distribuciones condicionales completas

(Ejercicio.) Las distribuciones condicionales completas para implementar un muestreador de Gibbs son:

  • La distribución condicional completa de \(\mu_j\) es \(\mu_j\mid\text{resto}\sim\textsf{N}(m_j,v_j^2)\) con: \[ m_j = \frac{\frac{1}{\tau^2}\,\boldsymbol{x}^{\textsf{T}}_j\boldsymbol{\beta} + \frac{n_j}{\sigma^2_j}\,\bar{y}_j}{\frac{1}{\tau^2} + \frac{n_j}{\sigma^2_j}} \qquad\text{y}\qquad v^2_j = \frac{1}{\frac{1}{\tau^2} + \frac{n_j}{\sigma^2_j}}\,. \] Alternativamente, la distribución condicional completa de \(\boldsymbol{\mu}\) es \(\boldsymbol{\mu}\mid\text{resto}\sim\textsf{N}(\boldsymbol{m}.\mathbf{V})\) con: \[ \boldsymbol{m} = \left(\tfrac{1}{\tau^2}\mathbf{I} + \mathbf{D}^{-1} \right)^{-1}\left(\tfrac{1}{\tau^2}\mathbf{X}\boldsymbol{\beta} + \mathbf{D}^{-1}\bar{\boldsymbol{y}}\right) \qquad\text{y}\qquad \mathbf{V} = \left(\tfrac{1}{\tau^2}\mathbf{I} + \mathbf{D}^{-1} \right)^{-1}\,. \] donde \(\mathbf{D}^{-1} = \textsf{diag}(n_1/\sigma^2_1,\ldots,n_m/\sigma^2_m)\).

  • La distribución condicional completa de \(\sigma^2_j\) es \(\sigma^2_j\mid\text{resto}\sim\textsf{GI}(a_j,b_j)\) con: \[ a_j = \frac{\nu+1}{2} \qquad\text{y}\qquad b_j = \frac{\nu\upsilon^2 + n_j(\bar{y}_j - \mu_j)^2}{2}\,. \]

  • La distribución condicional completa de \(\boldsymbol{\beta}\) es \(\boldsymbol{\beta}\mid\text{resto}\sim\textsf{N}(\boldsymbol{m},\mathbf{V})\) con: \[ \boldsymbol{m} = (\mathbf{X}^{\textsf{T}}\mathbf{X})^{-1}\mathbf{X}^{\textsf{T}}\boldsymbol{\mu} \qquad\text{y}\qquad \mathbf{V} = \tau^2(\mathbf{X}^{\textsf{T}}\mathbf{X})^{-1}\,. \]

  • La distribución condicional completa de \(\tau^2\) es \(\tau^2\mid\text{resto}\sim\textsf{GI}(a,b)\) con: \[ a = \frac{m}{2} \qquad\text{y}\qquad b = \frac{(\boldsymbol{\mu} - \mathbf{X}\boldsymbol{\beta})^{\textsf{T}}(\boldsymbol{\mu} - \mathbf{X}\boldsymbol{\beta})}{2}\,. \]

  • La distribución condicional completa de \(\nu\) es: \[ \log p(\nu\mid\text{resto}) \propto \frac{m\nu}{2}\log\left(\frac{\nu\upsilon^2}{2}\right) - m\log\Gamma\left(\frac{\nu}{2}\right) - \frac{\nu}{2}\sum_{j=1}^m\log\sigma^2_j - \nu\left( \lambda_0 + \frac{\upsilon^2}{2} \sum_{j=1}^m \frac{1}{\sigma^2_j} \right)\,. \]

  • La distribución condicional completa de \(\upsilon^2\) es \(\upsilon^2\mid\text{resto}\sim\textsf{G}(a,b)\) con: \[ a = \frac{\alpha_0 + m\nu}{2} \qquad\text{y}\qquad b = \frac{\beta_0 + \nu\sum_{j=1}^m\frac{1}{\sigma^2_j}}{2}\,. \]

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(\mathbf{\Theta}\mid\boldsymbol{y})\):

  1. Establecer el valor inicial de los parámetros.
  2. Actualizar el valor de los parámetros por medio de las distribuciones condicionales completas correspondientes.
  3. Almacenar el valor de los parámetros simulados.
sample_mu <- function(sj, X, nj, m, sig2, bet, tau2) 
{
  vj <- 1/(1/tau2 + nj/sig2)
  c(rnorm(n = m, mean = vj*(c(X%*%bet)/tau2 + sj/sig2), sd = sqrt(vj)))
}
sample_sig2 <- function(y, nj, m, mu, nu, ups2) 
{
  c(1/rgamma(n = m, shape = 0.5*(nu + 1), rate = 0.5*(nu*ups2 + nj*(y - mu)^2)))
}
sample_bet <- function(H, XtXi, mu, tau2)
{
  c(mvtnorm::rmvnorm(n = 1, mean = H%*%mu, sigma = tau2*XtXi))
}
sample_tau2 <- function(X, m, mu, bet)
{
  1/rgamma(n = 1, shape = 0.5*m, rate = 0.5*sum((mu - c(X%*%bet))^2))
}
sample_nu <- function(m, la0, sig2, nus, ups2) 
{
  lp <- 0.5*m*nus*log(0.5*nus*ups2) - 
        m*lgamma(0.5*nus) - 
        0.5*nus*sum(log(sig2)) - 
        nus*(la0 + 0.5*ups2*sum(1/sig2))
  sample(x = nus, size = 1, replace = F, prob = exp(lp - max(lp)))
}
sample_ups2 <- function(m, al0, be0, sig2, nu)
{
  rgamma(n = 1, shape = 0.5*(al0 + m*nu), rate = 0.5*(be0 + nu*sum(1/sig2)))
}
# MCMC
MCMC <- function (niter, nburn, nskip, y, X, nj, al0, be0, la0, nus, verbose = TRUE) 
{
  # Número iteraciones
  B <- nburn + nskip*niter
  ncat <- floor(0.1*B)
  
  # Constantes
  m    <- nrow(X)
  p    <- ncol(X)
  s2y  <- var(y)
  sj   <- c(nj*y)
  XtXi <- solve(t(X)%*%X)
  H    <- solve(t(X)%*%X)%*%t(X)
  
  # Inicialización
  mu   <- c(y)
  sig2 <- c(rep(s2y, m))
  bet  <- c(H%*%mu)
  tau2 <- c(sum((mu - X%*%bet)^2)/(m - p))
  nu   <- c(1)
  ups2 <- c(s2y)
  
  # Almacenamiento
  MU <- SIG2 <- BET <- TAU2 <- NU <- UPS2 <- LP <- NULL
  
  # Cadena
  set.seed(1234)
  for (b in 1:B) {
    # Actualizar
    mu   <- sample_mu  (sj, X, nj, m, sig2, bet, tau2)
    sig2 <- sample_sig2(y, nj, m, mu, nu, ups2)
    bet  <- sample_bet (H, XtXi, mu, tau2)
    tau2 <- sample_tau2(X, m, mu, bet)
    ups2 <- sample_ups2(m, al0, be0, sig2, nu)
    nu   <- sample_nu  (m, la0, sig2, nus, ups2)
    
    # Almacenar y log-verosimilitud
    if (b > nburn) {
      if (b%%nskip == 0) {
        MU   <- rbind(MU,   mu)
        SIG2 <- rbind(SIG2, sig2)
        BET  <- rbind(BET,  bet)
        TAU2 <- rbind(TAU2, tau2)
        NU   <- rbind(NU,   nu)
        UPS2 <- rbind(UPS2, ups2)
        LP   <- rbind(LP,   sum(dnorm(x = y, mean = mu, sd = sqrt(sig2/nj), log = T)))
      }
    
    }
    # progreso
    if (verbose) 
      if (b%%ncat == 0) 
        cat(100*round(b/B, 1), "% completado ... \n", sep = "")
  }
  
  # Retorno
  list(
       MU = MU, 
       SIG2 = SIG2, 
       BET = BET, 
       TAU2 = TAU2, 
       NU = NU, 
       UPS2 = UPS2, 
       LP = LP
  )
}

Selección de los hiperparámetros

Se escoge \(\lambda_0 = 0.5\) y \(\alpha_0 = 8\) y \(\beta_0 = 4.078888\) de forma que \[ \textsf{E}(\upsilon^2) = \frac{\alpha_0}{\beta_0} := \,s^2_{\bar{y}} \qquad\text{y}\qquad \textsf{CV}(\upsilon^2) = \sqrt{\frac{2}{\alpha_0}} := 0.5, \] donde \(s^2_{\bar{y}}=\frac{1}{m-1}\sum_{j=1}^m(\bar{y}_j-\bar{\bar{y}})^2\), con \(\bar{\bar{y}}=\frac{1}{m}\sum_{j=1}^m \bar{y}_j\), de manera que \(\sigma^2_j\) esté débilmente concentrado al rededor de \(s^2_{\bar{y}}\).

# Cálculo de los hiperparámetros de la distribución gamma inversa para upsilon^2
(al0 <- 2 / 0.5^2)
## [1] 8
(be0  <- al0 / var(y))
## [1] 4.078888
# Gráfico de la distribución previa para upsilon^2
par(mfrow = c(1, 1), mar = c(4.5, 3.5, 1.5, 1), mgp = c(2, 0.7, 0))

curve(expr = dgamma(x, shape = al0 / 2, rate = be0 / 2), from = 0, to = 10,
  lwd = 2, xlab = expression(upsilon^2), ylab = expression(p(upsilon^2)),
  main = expression("Distribución previa para " * upsilon^2))

# Hiperparámetro de la distribución previa
la0 <- 1

# Rango de posibles valores para nu
nus <- 1:20

# Cálculo de la densidad no normalizada
prior_nu <- exp(-la0 * nus)

# Gráfico de barras verticales de la distribución previa
par(mfrow = c(1, 1), mar = c(4.5, 3.5, 1.5, 1), mgp = c(2, 0.7, 0))

plot( x = nus, y = prior_nu, type = "h", lwd = 2,
  xlab = expression(nu), ylab = expression(p(nu)),
  main = "Distribución previa no normalizada para nu")

Matriz de diseño

Se define la matriz de diseño tal que \(\boldsymbol{x}_j=(x_{j,1},\ldots,x_{j,12})\) con

  • \(x_{j,1} = 1\).
  • \(x_{j,2} = 1\) si \(\text{tipo}_j = \text{eri}\) y \(x_{2,j} = 0\) en otro caso.
  • \(x_{j,3} = 1\) si \(\text{tipo}_j = \text{d.buoy}\) y \(x_{3,j} = 0\) en otro caso.
  • \(x_{j,4} = 1\) si \(\text{tipo}_j = \text{f.buoy}\) y \(x_{4,j} = 0\) en otro caso.
  • \(x_{j,5} = \text{latitud}_{j}\).
  • \(x_{j,6} = \text{latitud}_{j}\times x_{2_j}\).
  • \(x_{j.7} = \text{latitud}_{j}\times x_{3_j}\).
  • \(x_{j,8} = \text{latitud}_{j}\times x_{4_j}\).
  • \(x_{j,9} = \text{longitud}_{j}\).
  • \(x_{j,10} = \text{longitud}_{j}\times x_{2_j}\).
  • \(x_{j,11} = \text{longitud}_{j}\times x_{3_j}\).
  • \(x_{j,12} = \text{longitud}_{j}\times x_{4_j}\)

de forma que

  • Si \(\text{tipo}_j = \text{bucket}\), entonces: \[ \textsf{E}(\mu_j\mid\boldsymbol{x}_j,\boldsymbol{\beta}) = \beta_1 + \beta_5\,\text{latitud}_j + \beta_9\,\text{longitud}_j\,. \]
  • Si \(\text{tipo}_j = \text{eri}\), entonces: \[ \textsf{E}(\mu_j\mid\boldsymbol{x}_j,\boldsymbol{\beta}) = (\beta_1+\beta_2) + (\beta_5+\beta_6)\,\text{latitud}_j + (\beta_9+\beta_{10})\,\text{longitud}_j\,. \]
  • Si \(\text{tipo}_j = \text{d.buoy}\), entonces: \[ \textsf{E}(\mu_j\mid\boldsymbol{x}_j,\boldsymbol{\beta}) = (\beta_1+\beta_3) + (\beta_5+\beta_7)\,\text{latitud}_j + (\beta_9+\beta_{11})\,\text{longitud}_j\,. \]
  • Si \(\text{tipo}_j = \text{f.buoy}\), entonces: \[ \textsf{E}(\mu_j\mid\boldsymbol{x}_j,\boldsymbol{\beta}) = (\beta_1+\beta_4) + (\beta_5+\beta_8)\,\text{latitud}_j + (\beta_9+\beta_{12})\,\text{longitud}_j\,. \]
# Matriz de diseño
X5 <- as.matrix(df$lat)
X9 <- as.matrix(df$lon)
X2 <- as.matrix(c(rep(0,m1), rep(1,m2), rep(0,m3), rep(0,m4)))
X3 <- as.matrix(c(rep(0,m1), rep(0,m2), rep(1,m3), rep(0,m4)))
X4 <- as.matrix(c(rep(0,m1), rep(0,m2), rep(0,m3), rep(1,m4)))
X  <- cbind(1, X2, X3, X4, X5, X2*X5, X3*X5, X4*X5, X9, X2*X9, X3*X9, X4*X9)
m  <- nrow(X)
p  <- ncol(X)

# Dimensiones
(n <- nrow(X))
## [1] 86
(p <- ncol(X))
## [1] 12

Implementación

# Número de parámetros
2*m+p+3
## [1] 187
# Número de iteraciones
nskip <- 100
nburn <- 2000
niter <- 20000
# Semilla para reproducibilidad
set.seed(1234)

# Ajuste del modelo
muestras <- MCMC(niter, nburn, nskip, y, X, nj, al0, be0, la0, nus, verbose = T)

# Guardar muestras
save(muestras, file = "muestras_ejemplo_regresion_2.RData")

Convergencia

# Cargar muestras
load(file = "muestras_ejemplo_regresion_2.RData")

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     112   12238   17406   15390   20000   20452
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00011 0.00034 0.01344 0.01901 0.02024 0.24344

Inferencia: Medias específicas

Inferencia sobre \(\mu_j\):

# Función resumen: media, desviación estándar e IC al 95%
resumen_param <- function(x) {
  c(
    mean = mean(x),
    sd   = sd(x),
    q025 = quantile(x, 0.025),
    q975 = quantile(x, 0.975)
  )
}

# Número de muestras
B <- nrow(muestras$MU)

# Cálculo de ETA para cada iteración (CVs intra)
muestras$ETA <- matrix(NA, nrow = B, ncol = m)
for (b in 1:B) {
  muestras$ETA[b, ] <- sqrt(muestras$SIG2[b, ]/nj)/abs(muestras$MU[b, ])
}

# Estadísticas resumen para cada parámetro
mu_summ   <- apply(muestras$MU,   2, resumen_param)
sig2_summ <- apply(muestras$SIG2, 2, resumen_param)
beta_summ <- apply(muestras$BET,  2, resumen_param)
eta_summ  <- apply(muestras$ETA,  2, resumen_param)

# Codificación de color para los dispositivos
color <- c(
  rep(16, m1),  # gris oscuro: bucket
  rep(4,  m2),  # azul: eri
  rep(2,  m3),  # rojo: f.buoy
  rep(3,  m4)   # verde: d.buoy
)

Inferencia: Ratio de heterogeneidad entre vs dentro de grupos

Inferencia sobre \(\eta_j = \tau^2/(\sigma^2_j/n_j)\):

Inferencia: Componentes de varianza

Inferencia sobre \(\tau^2\) y \(\upsilon^2\):

\(\tau^2\) mide la variabilidad entre las medias de los dispositivos \(\mu_j\) respecto al predictor lineal \(\boldsymbol{x}_j^\top \boldsymbol{\beta}\). Cuanto menor sea \(\tau^2\), mayor será el encogimiento hacia el predictor lineal.

\(\upsilon^2\) es la escala promedio de las varianzas dentro de los dispositivos \(\sigma_j^2\). Controla el nivel general de dispersión dentro de los dispositivos.

Inferencia: Coeficientes de la regresión

Inferencia sobre \(\beta_j\):

El modelo indica diferencias moderadas entre tipos de dispositivos en relación con el intercepto, especialmente para f.buoy, cuyo diferencial respecto a bucket es más alto. Sin embargo, intervalos de credibilidad amplios muestran alta incertidumbre, por lo que no hay evidencia concluyente de efectos sistemáticos entre dispositivos.

Los efectos de latitud y longitud, así como sus diferenciales respecto al tipo de dispositivo, son en general pequeños e inciertos, con intervalos que incluyen cero.

Los resultados sugieren heterogeneidad entre grupos, pero sin un patrón claro asociado al tipo de dispositivo ni a la ubicación geográfica.

Encogimiento

El modelo aplica suavizamiento al estimar las medias grupales \(\mu_j\), haciendo que las estimaciones \(\hat{\mu}_j\) sean menos extremas que los promedios muestrales \(\bar{y}_j\). Este encogimiento es mayor en grupos con valores atípicos o mayor incertidumbre, y permite obtener estimaciones más estables y robustas, evitando el sobreajuste a la variabilidad muestral.

Estadísticos de prueba

# Semilla para reproducibilidad
set.seed(123)

# Función resumen de estadísticos
test_stats <- function(x) {
  c(
    mean     = mean(x),
    median   = median(x),
    sd       = sd(x),
    iqr      = diff(quantile(x, c(0.25, 0.75)))
  )
}

ts_display <- c("Media", "Mediana", "Desv. Estándar", "Rango Intercuartílico")

# Inicialización
TS <- NULL

# Generación de estadísticas simuladas
for (b in 1:niter) {
  mu   <- muestras$MU[b, ]
  sig2 <- muestras$SIG2[b, ]
  
  simulados <- rnorm(n = m, mean = mu, sd = sqrt(sig2 / nj))
  TS <- rbind(TS, test_stats(simulados))
}

# Estadísticos observados
TS0 <- test_stats(y)

Predicción

# Predicción en una grilla de latitud y longitud para cada tipo de dispositivo

# Semilla para reproducibilidad
set.seed(123)

# Selección de un grupo representativo (con mayor tamaño muestral) por tipo de dispositivo
g  <- df$type
j1 <- sample(which((g == "bucket")  & (nj == max(nj[g == "bucket"]))), 1)
j2 <- sample(which((g == "eri")     & (nj == max(nj[g == "eri"   ]))), 1)
j3 <- sample(which((g == "d.buoy")  & (nj == max(nj[g == "d.buoy"]))), 1)
j4 <-        which((g == "f.buoy")  & (nj == max(nj[g == "f.buoy"])))      # Único grupo

# Construcción de grilla espacial (latitud × longitud)
ng <- 25
gx <- seq(min(X[, 9]), max(X[, 9]), length.out = ng)  # Longitud
gy <- seq(min(X[, 5]), max(X[, 5]), length.out = ng)  # Latitud

# Inicialización de matrices de predicción posterior para cada tipo
z1 <- z2 <- z3 <- z4 <- array(NA, dim = c(ng, ng, 1000))

# Recorremos cada punto de la grilla
for (i in 1:ng) {
  for (j in 1:ng) {
    lon <- gx[i]
    lat <- gy[j]
    
    # Variables dummies para los 4 tipos de dispositivo
    type_dummy <- rbind(rep(0, 3), diag(3))  # bucket, eri, d.buoy, f.buoy

    # Matriz de diseño con interacciones latitud y longitud × tipo
    x <- cbind(1, type_dummy,          # intercepto y dummies
               lat, lat * type_dummy,  # latitud e interacción con tipo
               lon, lon * type_dummy)  # longitud e interacción con tipo

    # Muestreo predictivo posterior para cada tipo de dispositivo
    for (b in 1:1000) {
      bet  <- muestras$BET [b, ]
      sig2 <- muestras$SIG2[b, ]

      z1[i, j, b] <- rnorm(1, mean = sum(x[1, ] * bet), sd = sqrt(sig2[j1] / nj[j1]))
      z2[i, j, b] <- rnorm(1, mean = sum(x[2, ] * bet), sd = sqrt(sig2[j2] / nj[j2]))
      z3[i, j, b] <- rnorm(1, mean = sum(x[3, ] * bet), sd = sqrt(sig2[j3] / nj[j3]))
      z4[i, j, b] <- rnorm(1, mean = sum(x[4, ] * bet), sd = sqrt(sig2[j4] / nj[j4]))
    }
  }
}
filled.contour(x = gx, 
  y = gy, 
  z = apply(z1, c(1, 2), mean),
  xlab = "Latitud", 
  ylab = "Longitud", 
  main = "bucket",
  las  = 0)

# Estimación de la media posterior por punto de grilla para cada tipo de dispositivo
zz1 <- as.vector(apply(z1, c(1, 2), mean))
zz2 <- as.vector(apply(z2, c(1, 2), mean))
zz3 <- as.vector(apply(z3, c(1, 2), mean))
zz4 <- as.vector(apply(z4, c(1, 2), mean))
# Temperatura promedio global (ponderada por tamaño del grupo)
yb <- sum(y * nj) / sum(nj)

# Tabla de resumen: media predicha, RMSE, número de dispositivos y total de observaciones por tipo
tab <- rbind(
  c(mean(zz1), sqrt(mean((zz1 - yb)^2)), m1, sum(df1$n)),
  c(mean(zz2), sqrt(mean((zz2 - yb)^2)), m2, sum(df2$n)),
  c(mean(zz3), sqrt(mean((zz3 - yb)^2)), m3, sum(df3$n)),
  c(mean(zz4), sqrt(mean((zz4 - yb)^2)), m4, sum(df4$n))
)

# Asignar nombres a filas y columnas
colnames(tab) <- c("Media", "RMSE", "n disp", "n total")
rownames(tab) <- c("bucket", "eri", "d.buoy", "f.buoy")

# Mostrar tabla redondeada
round(tab, 2)
##        Media RMSE n disp n total
## bucket 19.40 0.57     36      36
## eri    19.74 0.53     35      39
## d.buoy 19.93 0.36     10      21
## f.buoy 22.66 3.25      5     240
# Comparación con valor de referencia de la NOAA para 2012
# Temperatura promedio según NOAA
yb <- 22.44

# Tabla: media predicha, RMSE respecto a NOAA, número de dispositivos, total de observaciones
tab <- rbind(
  c(mean(zz1), sqrt(mean((zz1 - yb)^2)), m1, sum(df1$n)),
  c(mean(zz2), sqrt(mean((zz2 - yb)^2)), m2, sum(df2$n)),
  c(mean(zz3), sqrt(mean((zz3 - yb)^2)), m3, sum(df3$n)),
  c(mean(zz4), sqrt(mean((zz4 - yb)^2)), m4, sum(df4$n))
)

# Asignar nombres de columnas y filas
colnames(tab) <- c("Media", "RMSE", "n disp", "n total")
rownames(tab) <- c("bucket", "eri", "d.buoy", "f.buoy")

# Mostrar tabla redondeada
round(tab, 2)
##        Media RMSE n disp n total
## bucket 19.40 3.04     36      36
## eri    19.74 2.74     35      39
## d.buoy 19.93 2.54     10      21
## f.buoy 22.66 1.83      5     240

Ejercicios

  • Cargue el conjunto de datos Boston Housing en R que se encuentra disponible en la libreria MASS. La variable de respuesta es medv, que corresponde al valor medio de las viviendas ocupadas por sus propietarios (en miles de dólares). Las otras 13 variables son covariables que describen las características del vecindario.

    1. Ajuste un modelo de regresión lineal con una previa difusa independiente para los coeficientes de regresión. Verifique que el muestreador MCMC haya convergido y resuma la distribución posterior de todos los coeficientes de la regresión.
    2. Reajuste el modelo utilizando previas Ridge y Lasso para los coeficientes de regresión y compare los resultados respecto al análisis con la previa Normal no informativa.
    3. Ajuste los modelos de a. y b. utilizando únicamente las primeras 500 observaciones y calcule la distribución predictiva posterior para las últimas 6 observaciones. Grafique la distribución predictiva posterior frente al valor real para estas 6 observaciones y comente si las predicciones son razonables.
  • Considere la base de datos de diabetes descrita en la Sección 9.3 de Hoff (2009), que incluye 10 medidas basales \(x_1, \ldots, x_{10}\) para un grupo de 442 pacientes diabéticos, junto con una medida de progresión de la enfermedad \(y\) registrada un año después. Los datos pueden descargarse desde este enlace en los archivos yX.diabetes.train y yX.diabetes.test.

    El objetivo es construir un modelo predictivo para \(y\) a partir de \(x_1, \ldots, x_{10}\), con todas las variables estandarizadas. Aunque un modelo de regresión con diez variables no sería excesivamente complejo, se sospecha que la relación entre \(y\) y las \(x_j\) podría no ser lineal. Por ello, se sugiere incluir términos de segundo orden, como \(x_j^2\) y \(x_j x_k\), para mejorar la capacidad predictiva.

    En consecuencia, las variables regresoras comprenden diez efectos principales \(x_j\), \(\binom{10}{2} = 45\) interacciones \(x_j x_k\) y nueve términos cuadráticos \(x_j^2\) (no se considera \(x_2^2\) porque \(x_2 = \text{sexo}\) es binaria y cumple \(x_2 = x_2^2\)). Esto resulta en un total de \(p = 64\) variables regresoras, sin incluir intercepto debido a la estandarización de todas las variables.

    Se considera el modelo de regresión \(\boldsymbol{y} \mid \mathbf{X}, \boldsymbol{\beta}, \sigma^2 \sim \textsf{N}_n(\mathbf{X}\boldsymbol{\beta}, \sigma^2\mathbf{I})\), donde \(\boldsymbol{y}\) es un vector \(n \times 1\) con los valores de la variable respuesta, \(\mathbf{X}\) es una matriz \(n \times p\) con las variables regresoras, y \(\boldsymbol{\beta}\) es un vector \(p \times 1\) de parámetros desconocidos.

    Para evaluar el modelo, los 442 individuos con diabetes se dividieron aleatoriamente en 342 para entrenamiento y 100 para prueba, obteniendo así los conjuntos \((\boldsymbol{y}_{\text{train}}, \mathbf{X}_{\text{train}})\) y \((\boldsymbol{y}_{\text{test}}, \mathbf{X}_{\text{test}})\). El ajuste se realiza con el conjunto de entrenamiento y, a partir de los coeficientes estimados \(\hat{\boldsymbol{\beta}} = \textsf{E}(\boldsymbol{\beta} \mid \boldsymbol{y}_{\text{train}})\), se generan las predicciones para el conjunto de prueba \(\hat{\boldsymbol{y}}_{\text{test}} = \mathbf{X}_{\text{test}} \hat{\boldsymbol{\beta}}\). El rendimiento predictivo se determina comparando \(\hat{\boldsymbol{y}}_{\text{test}}\) con \(\boldsymbol{y}_{\text{test}}\) mediante una métrica adecuada.

    Considere los modelos de regresión con previa unitaria, regresión con previa \(g\), regresión con previa difusa independiente, regresión Ridge, y regresión Lasso.

    1. Ajuste cada modelo utilizando los datos de entrenamiento \((\boldsymbol{y}_{\text{train}},\mathbf{X}_{\text{train}})\).
    2. Para cada modelo, calcule \(\hat{\boldsymbol{y}}_{\text{test}} = \mathbf{X}_{\text{test}}\hat{\boldsymbol{\beta}}\) usando los coeficientes de regresión estimados \(\hat{\boldsymbol{\beta}} = \textsf{E}(\boldsymbol{\beta}\mid\boldsymbol{y}_{\text{train}})\). Grafique \(\hat{y}_{\text{test}}\) frente \(y_{\text{test}}\) y calcule el el coeficiente de determinación, la raíz del error cuadrático medio, y el error absoluto medio.
    3. Para cada modelo, establezca la bondad de ajuste usando la media y la desviación estándar como estadísticos de pruebas.
    4. Para cada modelo, calcule el DIC y el WAIC.
    5. Interprete los resultados obtenidos en los numerales anteriores.

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.