13.6 Practice

The practice problems for Chapter 13 ask you to fit varying slopes models, interpret covariance matrices, and explore Gaussian process priors.

Practice

Easy problems

13E1. What is a varying slope? How does it differ from a varying intercept?

13E2. What does the correlation $\rho$ between intercepts and slopes tell us?

13E3. What is the LKJ prior? What happens when $\eta$ increases?

13E4. What is a Gaussian process? What does the length scale $\ell$ control?

Medium problems

13M1. Fit a varying slopes model. Interpret the covariance matrix and the correlation.

# 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")

Click Run to execute this code in your browser with WebR

13M2. Simulate Gaussian process data and fit a GP regression. Explore the effect of the length scale.

# 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")

Click Run to execute this code in your browser with WebR

Hard problems

13H1. Fit a varying slopes model with a negative correlation between intercepts and slopes. Verify that the model recovers the correlation.

# 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")

Click Run to execute this code in your browser with WebR

13H2. Fit a Gaussian process to spatial data. Use different kernels and compare their predictions.

13H3. Compare a varying slopes model with and without the correlation parameter. Does allowing $\rho \neq 0$ improve the model? Use AIC or WAIC.

Summary

Measurement error