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}\,. \]
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*} \]
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).
Motivación:
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)\).
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:
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:
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!):
Este método aprovecha la conjugación inherente a la previa \(g\) para facilitar simulaciones directas de la distribución posterior.
Referencias:
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\).
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:
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:
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:
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
## [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
# 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
## beta1 beta2 beta3 beta4
## 0.001 0.768 0.999 0.330
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
}
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.
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.
(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}\,. \]
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.
(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\,. \]
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.
# 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
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
))
}
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
))
}
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
))
}
# 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%)
# 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")
# 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$")
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)
}
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.")
RMSE | MAE | |
---|---|---|
OLS | 0.458 | 0.377 |
g-prior | 0.449 | 0.371 |
Ridge | 0.341 | 0.289 |
Lasso | 0.427 | 0.353 |
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.
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\).
## [1] 86 6
## 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}}\).
## [1] 86
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.00 18.70 19.65 19.56 20.60 22.00
## nj
## 1 2 3 23 24 37 62 94
## 70 7 4 1 1 1 1 1
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\).
(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*} \]
(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}\,. \]
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})\):
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
)
}
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}}\).
## [1] 8
## [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")
Se define la matriz de diseño tal que \(\boldsymbol{x}_j=(x_{j,1},\ldots,x_{j,12})\) con
de forma que
# 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
## [1] 12
## [1] 187
## 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 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 sobre \(\eta_j = \tau^2/(\sigma^2_j/n_j)\):
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 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.
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.
# 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 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
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.
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.
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.