13.6 Práctica
Los problemas de práctica del Capítulo 13 le piden ajustar modelos de pendientes variables, interpretar matrices de covarianza y explorar distribuciones a priori de procesos Gaussianos.
Problemas fáciles
13E1. ¿Qué es una pendiente variable? ¿En qué se diferencia de un intercepto variable?
13E2. ¿Qué nos dice la correlación entre interceptos y pendientes?
13E3. ¿Qué es la distribución a priori LKJ? ¿Qué pasa cuando aumenta?
13E4. ¿Qué es un proceso Gaussiano? ¿Qué controla la escala de longitud ?
Problemas medios
13M1. Ajuste un modelo de pendientes variables. Interprete la matriz de covarianza y la correlación.
# 13M1: Varying slopes model
set.seed(42)
n_groups <- 15
n_per <- 20
a_bar <- 2
b_bar <- 0.5
sigma_a <- 1.0
sigma_b <- 0.3
rho <- -0.4
library(MASS)
Sigma <- matrix(c(sigma_a^2, rho*sigma_a*sigma_b,
rho*sigma_a*sigma_b, sigma_b^2), nrow = 2)
params <- mvrnorm(n_groups, c(0, 0), Sigma)
d <- do.call(rbind, lapply(1:n_groups, function(j) {
x <- runif(n_per, 0, 10)
y <- a_bar + params[j, 1] + (b_bar + params[j, 2]) * x + rnorm(n_per, 0, 0.5)
data.frame(group = factor(j), x = x, y = y)
}))
library(lme4)
m <- lmer(y ~ x + (1 + x | group), data = d)
cat("Varying slopes model:\n")
cat(" Fixed: intercept =", round(fixef(m)[1], 3),
" slope =", round(fixef(m)[2], 3), "\n")
cat(" SD: intercept =", round(attr(VarCorr(m)$group, "stddev")[1], 3),
" slope =", round(attr(VarCorr(m)$group, "stddev")[2], 3), "\n")
cat(" Correlation:", round(attr(VarCorr(m)$group, "correlation")[1, 2], 3), "\n")
cat("\n True correlation was", rho, "\n")13M2. Simule datos de proceso Gaussiano y ajuste una regresión con PG. Explore el efecto de la escala de longitud.
# 13M2: Gaussian process regression
set.seed(42)
x <- seq(0, 10, length.out = 40)
f_true <- 2 * sin(x) + 0.3 * x
y <- f_true + rnorm(length(x), 0, 0.3)
# Try different length scales
for (ell in c(0.5, 1.5, 3)) {
sigma*f <- 2
sigma_n <- 0.3
K <- sigma_f^2 * exp(-outer(x, x, "-")^2 / (2 _ ell^2))
K_noise <- K + sigma_n^2 * diag(length(x))
L <- chol(K_noise)
alpha <- backsolve(L, forwardsolve(t(L), y))
x*new <- seq(0, 10, length.out = 100)
K_star <- sigma_f^2 * exp(-outer(x*new, x, "-")^2 / (2 * ell^2))
f_pred <- K_star %*% alpha
rmse <- sqrt(mean((f*pred - (2 * sin(x*new) + 0.3 * x_new))^2))
cat(sprintf(" ell = %.1f: RMSE = %.3f\n", ell, rmse))
}
cat("\n Small ell: wiggly fit (overfits noise)\n")
cat(" Large ell: too smooth (underfits)\n")
cat(" Best ell: balances fit and smoothness\n")Problemas difíciles
13H1. Ajuste un modelo de pendientes variables con una correlación negativa entre interceptos y pendientes. Verifique que el modelo recupera la correlación.
# 13H1: Recover correlation from varying slopes model
set.seed(42)
n_groups <- 20
n_per <- 25
a_bar <- 3
b_bar <- 0.4
sigma_a <- 1.2
sigma_b <- 0.4
rho_true <- -0.7
library(MASS)
Sigma <- matrix(c(sigma_a^2, rho_true*sigma_a*sigma_b,
rho_true*sigma_a*sigma_b, sigma_b^2), nrow = 2)
params <- mvrnorm(n_groups, c(0, 0), Sigma)
d <- do.call(rbind, lapply(1:n_groups, function(j) {
x <- runif(n_per, 0, 10)
y <- a_bar + params[j, 1] + (b_bar + params[j, 2]) * x + rnorm(n_per, 0, 0.5)
data.frame(group = factor(j), x = x, y = y)
}))
library(lme4)
m <- lmer(y ~ x + (1 + x | group), data = d)
rho_est <- attr(VarCorr(m)$group, "correlation")[1, 2]
cat("Correlation recovery:\n")
cat(" True rho:", rho_true, "\n")
cat(" Estimated rho:", round(rho_est, 3), "\n")
cat(" SD intercept: true =", sigma_a, " est =",
round(attr(VarCorr(m)$group, "stddev")[1], 3), "\n")
cat(" SD slope: true =", sigma_b, " est =",
round(attr(VarCorr(m)$group, "stddev")[2], 3), "\n")13H2. Ajuste un proceso Gaussiano a datos espaciales. Use diferentes kernels y compare sus predicciones.
13H3. Compare un modelo de pendientes variables con y sin el parámetro de correlación. ¿Mejora el modelo permitir ? Use AIC o WAIC.