Modelo

Este modelo asume que la variable de respuesta \(\boldsymbol{y}\) sigue una distribución normal multivariante, condicionada a un conjunto de covariables \(\mathbf{X}\) y a parámetros desconocidos \(\boldsymbol{\beta}\), \(\sigma^2\) y \(\rho\).

La dependencia entre las observaciones se representa mediante una matriz de covarianza \(\mathbf{C}_\rho\), cuya estructura se basa en un proceso autorregresivo de primer orden (AR(1)).

Esta especificación captura patrones de correlación espacial o temporal, donde la fuerza de la asociación disminuye exponencialmente a medida que aumenta la distancia entre observaciones. El modelo resulta especialmente adecuado para analizar datos secuenciales o espacialmente organizados.

Verosimilitud

La distribución muestral está dada por:
\[ \boldsymbol{y}\mid\mathbf{X},\boldsymbol{\beta},\sigma^2,\rho \sim \textsf{N}\left(\mathbf{X}\boldsymbol{\beta}, \sigma^2 \mathbf{C}_\rho \right), \]
donde \(\mathbf{C}_\rho\) representa una matriz de covarianza con estructura autorregresiva de primer orden:
\[ \mathbf{C}_\rho = \begin{bmatrix} 1 & \rho & \rho^{2} & \cdots & \rho^{n-1} \\ \rho & 1 & \rho & \cdots & \rho^{n-2} \\ \rho^{2} & \rho & 1 & \cdots & \rho^{n-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho^{n-1} & \rho^{n-2} & \rho^{n-3} & \cdots & 1 \end{bmatrix}. \]
Este tipo de matriz garantiza que las correlaciones entre observaciones adyacentes sean más fuertes y disminuyan progresivamente a medida que aumenta la distancia entre ellas.

Distribución previa

Se asumen las siguientes distribuciones previas para los parámetros:
\[ p(\boldsymbol{\beta}, \sigma^2, \rho) = p(\boldsymbol{\beta}) \, p(\sigma^2) \, p(\rho), \]
con las siguientes especificaciones:

  • Para los coeficientes \(\boldsymbol{\beta}\):
    \[ \boldsymbol{\beta} \sim \textsf{N}(\boldsymbol{\beta}_0, \boldsymbol{\Sigma}_0). \]
  • Para la varianza residual \(\sigma^2\):
    \[ \sigma^2 \sim \textsf{GI}\left(\frac{\nu_0}{2}, \frac{\nu_0 \, \sigma^2_0}{2}\right). \]
  • Para el parámetro de autocorrelación \(\rho\):
    \[ \rho \sim \textsf{U}(a, b), \]

Parámetros e hiperparámetros

  • Parámetros del modelo:
    \[ \boldsymbol{\theta} = (\boldsymbol{\beta}, \sigma^2, \rho). \]

  • Hiperparámetros:
    \[ (\boldsymbol{\beta}_0, \boldsymbol{\Sigma}_0, \nu_0, \sigma^2_0, a, b). \]

Estimación

e propone construir un muestreador de Gibbs para generar muestras de la distribución posterior \(p(\boldsymbol{\theta} \mid \boldsymbol{y})\), permitiendo la estimación bayesiana de los parámetros del modelo.

Distribuciones condicionales completas

  1. Para los coeficientes \(\boldsymbol{\beta}\):
    \[ \boldsymbol{\beta} \mid \text{resto} \sim \textsf{N}\left(\boldsymbol{\beta}_{n}, \boldsymbol{\Sigma}_{n}\right), \]
    donde:
    \[ \boldsymbol{\Sigma}_{n} = \left(\boldsymbol{\Sigma}_{0}^{-1} + \frac{1}{\sigma^{2}} \mathbf{X}^{T} \mathbf{C}_{\rho}^{-1} \mathbf{X}\right)^{-1}, \]
    \[ \boldsymbol{\beta}_{n} = \boldsymbol{\Sigma}_{n}\left(\boldsymbol{\Sigma}_{0}^{-1} \boldsymbol{\beta}_{0} + \frac{1}{\sigma^{2}} \mathbf{X}^{T} \mathbf{C}_{\rho}^{-1} \boldsymbol{y}\right). \]

  2. Para la varianza residual \(\sigma^2\):
    \[ \sigma^{2} \mid \text{resto} \sim \textsf{IG}\left(\frac{\nu_{0} + n}{2}, \frac{\nu_{0} \, \sigma_{0}^{2} + \mathrm{SSR}_{\rho}}{2}\right), \]
    donde:
    \[ \mathrm{SSR}_{\rho} = (\boldsymbol{y} - \mathbf{X}\boldsymbol{\beta})^{\textsf{T}} \mathbf{C}_{\rho}^{-1} (\boldsymbol{y} - \mathbf{X}\boldsymbol{\beta}). \]

  3. La distribución condicional completa (DCC) del parámetro de autocorrelación \(\rho\) no posee una expresión analítica cerrada.

Algoritmo de Metropolis para muestrear \(\rho\)

Se utiliza un muestreador de Metropolis en un espacio transformado para garantizar que \(\rho \in (a, b)\). Para ello, se define una transformación que mapea \(\rho\) al espacio real \(\mathbb{R}\), lo que facilita el muestreo.

Transformación directa:
Mapear \(\rho \in (a, b)\) a \(\gamma \in \mathbb{R}\):
\[ \gamma = \log\left( \frac{\rho - a}{b - \rho} \right). \]

Transformación inversa:
Recuperar \(\rho\) desde \(\gamma\):
\[ \rho = \frac{a + b \exp(\gamma)}{1 + \exp(\gamma)}. \]

Jacobiano de la transformación:
El término de ajuste requerido para la densidad posterior es:
\[ \frac{d\rho}{d\gamma} = \frac{(b - a) \exp(\gamma)}{(1 + \exp(\gamma))^2}. \]

Algoritmo de Metropolis

  1. Proponer un valor candidato:
    Simular una propuesta \(\gamma^*\) a partir de una distribución normal centrada en el estado actual:
    \[ \gamma^* \sim \textsf{N}(\gamma^{(t-1)}, \delta^2), \]
    donde \(\delta^2\) es un parámetro de ajuste (tuning parameter) que controla la varianza de la propuesta.

  2. Transformar al espacio original:
    Calcular el valor propuesto para \(\rho^*\):
    \[ \rho^* = \frac{a + b \exp(\gamma^*)}{1 + \exp(\gamma^*)}. \]

  3. Calcular la probabilidad de aceptación:
    La probabilidad de aceptación incorpora la función objetivo y el Jacobiano:
    \[ r = \exp\left[ \log p(\rho^* \mid \text{resto}) - \log p(\rho^{(t-1)} \mid \text{resto}) + \log \left( \frac{d\rho^*}{d\gamma^*} \right) - \log \left( \frac{d\rho^{(t-1)}}{d\gamma^{(t-1)}} \right) \right]. \]

  4. Aceptar o rechazar la propuesta:
    Aceptar con probabilidad \(\min(1, r)\):
    \[ \gamma^{(t)} = \begin{cases} \gamma^* & \text{con probabilidad } \min(1, r), \\ \gamma^{(t-1)} & \text{en caso contrario.} \end{cases} \]

  5. Actualizar \(\rho^{(t)}\):
    Transformar el nuevo valor aceptado de \(\gamma^{(t)}\) al intervalo \((a, b)\):
    \[ \rho^{(t)} = \frac{a + b \exp(\gamma^{(t)})}{1 + \exp(\gamma^{(t)})}. \]

Hiperparámetros

Media previa \(\boldsymbol{\beta}_0\)

  • Elección estándar: Usar \(\boldsymbol{\beta}_0 = \mathbf{0}\) para priorizar un modelo centrado en valores cercanos a cero, lo que introduce un sesgo mínimo en ausencia de información previa.

Covarianza previa \(\boldsymbol{\Sigma}_0\)

  • No informativa: Utilizar una matriz diagonal con valores grandes para reflejar alta incertidumbre y permitir que los datos dominen el ajuste.

Parámetros para la varianza residual \(\nu_0, \sigma^2_0\)

  • \(\nu_0\) (grados de libertad): Usar valores pequeños (\(\nu_0 = 1\) o \(2\)) para imponer una distribución poco informativa, permitiendo que los datos controlen la estimación.
  • \(\sigma^2_0\) (escala inicial): Estimar la varianza empírica de los residuos a partir de un modelo de mínimos cuadrados ordinarios (OLS) como referencia:
    \[ \sigma^2_0 = \frac{1}{n - p}(y - X \hat{\beta})^T(y - X \hat{\beta}), \]
    donde \(p\) representa el número de parámetros. Este enfoque garantiza una especificación inicial coherente con la estructura observada en los datos.

Parámetros para el parámetro de autocorrelación (\(a, b\))

  • Correlación positiva esperada: Fijar \(a = 0\) y \(b = 1\).
  • Correlaciones negativas permitidas: Usar \(a = -1\) y \(b = 1\) para abarcar relaciones inversas potenciales.

Parámetro inicial y ajuste adaptativo \(\delta^2\)

  • Inicialización: Establecer \(\delta^2 = 0.1\) o \(0.01\) para facilitar propuestas razonables en las primeras iteraciones.
  • Ajuste adaptativo: Aprovechar el método exponencial ya implementado para ajustar dinámicamente \(\delta^2\) y mantener una tasa de aceptación en el rango deseado (30%–50%).

Regresión lineal con errores correlacionados

Los análisis de núcleos de hielo extraídos de la Antártida han permitido a los científicos reconstruir las condiciones climáticas y atmosféricas a lo largo de cientos de miles de años.

Estos registros, obtenidos mediante perforaciones profundas en capas de hielo acumuladas durante milenios, preservan burbujas de aire atrapadas que proporcionan información directa sobre la composición atmosférica pasada, incluyendo niveles de dióxido de carbono y metano.

Estudios como el realizado en el núcleo de hielo de Vostok han revelado patrones de variabilidad climática, ciclos glaciales e interglaciales, y la relación entre los gases de efecto invernadero y la temperatura global.

Petit, J. R., et al. (1999). Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica. Nature, 399(6735), 429–436.

El conjunto de datos consta de 200 observaciones de temperatura registradas en intervalos de tiempo aproximadamente equidistantes, con un lapso de alrededor de 2,000 años entre mediciones consecutivas.

La temperatura se expresa como la diferencia relativa respecto a la temperatura actual, medida en grados Celsius, mientras que la concentración de CO\(_2\) (dióxido de carbono) se reporta en partes por millón (ppm).

El objetivo es modelar la temperatura como una función de la concentración de CO\(_2\), explorando la relación entre ambas variables para inferir patrones climáticos históricos y evaluar posibles vínculos entre los niveles de gases de efecto invernadero y el cambio de temperatura.

Data

  • \(y_i\,\,\): CO2 registrada en el año \(i\), para \(i=1,\ldots,n\).
  • \(x_{i}\,\): temperatura registrada en el año \(i\).
# datos
load("dct.RData")
head(dct, n = 5)
##      year   co2  tmp
## 1 -419095 277.6 0.94
## 2 -416872 286.9 1.36
## 3 -414692 285.5 1.07
## 4 -413062 285.5 1.08
## 5 -410979 281.2 1.27
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# ajuste del modelo clásico
lmfit <- lm(dct$tmp ~ dct$co2)
summary(lmfit)
## 
## Call:
## lm(formula = dct$tmp ~ dct$co2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2351 -0.9715  0.0082  1.1544  4.4754 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -23.024140   0.879543  -26.18   <2e-16 ***
## dct$co2       0.079853   0.003834   20.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.533 on 198 degrees of freedom
## Multiple R-squared:  0.6866, Adjusted R-squared:  0.6851 
## F-statistic: 433.8 on 1 and 198 DF,  p-value: < 2.2e-16
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Ajuste del modelo

# datos
n  <- dim(dct)[1]
y  <- dct[,3]
X  <- cbind(rep(1,n), dct[,2])
# Algoritmo de Gibbs para el modelo con estructura AR(1)

# Librerías necesarias
suppressMessages(suppressWarnings(library(nlme)))
suppressMessages(suppressWarnings(library(MASS)))  # Para mvrnorm

# Datos y configuraciones iniciales
DY <- abs(outer(1:n, 1:n, "-"))  # Matriz de distancias
lmfit <- lm(y ~ -1 + X)
fit.gls <- gls(y ~ X[, 2], correlation = corARMA(p = 1), method = "ML")

# Valores iniciales
beta <- lmfit$coef
s2 <- summary(lmfit)$sigma^2
rho <- acf(lmfit$res, plot = FALSE)$acf[2]

# Hiperparámetros
nu0 <- 1
s20 <- summary(lmfit)$sigma^2
iSig0 <- diag(1 / 1000, nrow = 2)

# Configuración del algoritmo
n_iter <- 55000
n_burn <- 5000
n_thin <- 10
n_cat <- floor(n_iter / 10)  # Mostrar progreso cada 10%
samples <- NULL
ac <- 0

# Parámetro adaptativo
delta2 <- 0.1
accept_rate_target <- 0.44
adapt_rate <- 0.05

# Algoritmo de Gibbs
set.seed(123)
for (b in 1:n_iter) {
  # Simular beta
  Cor <- rho^DY
  iCor <- solve(Cor)
  V_beta <- solve(t(X) %*% iCor %*% X / s2 + iSig0)
  E_beta <- V_beta %*% (t(X) %*% iCor %*% y / s2)
  beta <- c(mvrnorm(1, E_beta, V_beta))  # Usar mvrnorm de MASS

  # Simular sigma^2
  s2 <- 1 / rgamma(1, (nu0 + n) / 2, (nu0 * s20 + t(y - X %*% beta) %*% iCor %*% (y - X %*% beta)) / 2)

  # Simular rho (Metropolis adaptativo)
  # Paso 1: Propuesta en espacio transformado
  gamma_c <- log(rho / (1 - rho))  # Transformar a espacio real
  gamma_p <- rnorm(1, gamma_c, sqrt(delta2))  # Propuesta normal adaptativa
  rho_p <- exp(gamma_p) / (1 + exp(gamma_p))  # Transformar de regreso

  # Paso 2: Probabilidad de aceptación con Jacobiano
  log_jacobian_c <- log(rho) + log(1 - rho)
  log_jacobian_p <- log(rho_p) + log(1 - rho_p)

  r <- exp(-0.5 * (det(rho_p^DY) - det(rho^DY) +
    sum(diag((y - X %*% beta) %*% t(y - X %*% beta) %*% (solve(rho_p^DY) - solve(rho^DY)))) / s2) +
    (log_jacobian_p - log_jacobian_c))

  # Paso 3: Aceptar o rechazar la propuesta
  if (runif(1) < r) {
    rho <- rho_p
    ac <- ac + 1
  }

  # Adaptación de delta2
  accept_rate <- ac / b
  delta2 <- delta2 * exp(adapt_rate * (accept_rate - accept_rate_target))

  # Almacenar resultados después de burn-in y thinning
  if (b > n_burn && (b - n_burn) %% n_thin == 0) {
    samples <- rbind(samples, c(beta, s2, rho))
  }

  # Mostrar progreso
  if (b %% n_cat == 0) {
    cat(sprintf("Iteración: %d, Tasa de aceptación: %.2f%%, delta2: %.3f, beta[0]: %.3f, beta[1]: %.3f, sigma^2: %.3f, rho: %.3f\n",
                b, (ac / b) * 100, delta2, beta[1], beta[2], s2, rho))
  }
}

# Resultados
colnames(samples) <- c("beta[0]", "beta[1]", "sigma^2", "rho")
samples <- as.data.frame(samples)
save(samples, file = "samples_dtc.RData")

Convergencia

# Calcular tamaños de muestra efectiva
coda::effectiveSize(samples)
##  beta[0]  beta[1]  sigma^2      rho 
## 5000.000 5000.000 5000.000 5055.743

Inferencia

# Calcular estadísticas para cada parámetro
stats <- data.frame(
  Parámetro = colnames(samples),
  Estimación = round(apply(samples, 2, mean), 3),
  CV = round(abs(apply(samples, 2, sd) / apply(samples, 2, mean)) * 100, 3),
  Límite_Inferior = round(apply(samples, 2, function(x) quantile(x, 0.025)), 3),
  Límite_Superior = round(apply(samples, 2, function(x) quantile(x, 0.975)), 3)
)

# Mostrar la tabla
print(stats, row.names = FALSE)
##  Parámetro Estimación     CV Límite_Inferior Límite_Superior
##     beta_0    -21.458  5.365         -23.655         -19.171
##     beta_1      0.073  6.869           0.063           0.083
##     sigma2      2.029 10.161           1.660           2.469
##        rho      0.304 16.165           0.207           0.399

  • \(\beta_0 = -21.458\): Representa la temperatura promedio cuando la concentración de CO\(_2\) es cero. Aunque este valor no es físicamente realista, sirve como punto de referencia para evaluar cambios en la temperatura conforme aumentan los niveles de CO\(_2\).

  • \(\beta_1 = 0.073\): Indica que un incremento de 1 ppm en la concentración de CO\(_2\) está asociado con un aumento promedio de 0.073 °C en la temperatura. Esto evidencia una relación positiva entre el CO\(_2\) y la temperatura, apoyando la influencia de los gases de efecto invernadero en el calentamiento global.

  • \(\sigma^2 = 2.029\): Representa la variabilidad no explicada en la temperatura. Muestra que hay factores adicionales no incluidos en el modelo que influyen en las mediciones, sugiriendo que el CO\(_2\) explica parte, pero no toda, la variabilidad observada.

  • \(\rho = 0.304\): Refleja la dependencia temporal entre observaciones separadas por 2,000 años. Indica que las temperaturas consecutivas tienden a estar moderadamente correlacionadas, lo cual es consistente con la persistencia climática observada en los registros históricos.

Referencias