13.1 Varying slopes by construction

Varying slopes extend varying intercepts by allowing the effect of a predictor to differ across groups. The model estimates a covariance matrix of intercepts and slopes.

13.1.1 The varying slopes model

Instead of only varying intercepts, we now let both the intercept and slope vary by group:

$\begin{pmatrix} \gamma_j \\ \beta_j \end{pmatrix} \sim \text{MVNormal}\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_\gamma^2 & \rho\sigma_\gamma\sigma_\beta \\ \rho\sigma_\gamma\sigma_\beta & \sigma_\beta^2 \end{pmatrix}\right)$

The key addition is the correlation $\rho$ between intercepts and slopes — groups that start high may also change faster or slower.

# Varying slopes: cafe example
# Each cafe has its own intercept (morning visits) and slope (afternoon effect)

set.seed(42)
n_cafes <- 20

# Population parameters

a_bar <- 3.5 # average morning visits
b_bar <- -0.9 # average afternoon effect (fewer visits)
sigma_a <- 1.0 # between-cafe SD for intercepts
sigma_b <- 0.5 # between-cafe SD for slopes
rho <- -0.6 # correlation: cafes with more morning visits lose more in afternoon

# Draw cafe-specific parameters from MVN

library(MASS)
Sigma <- matrix(c(sigma*a^2, rho * sigma*a * sigma*b,
rho * sigma*a * sigma_b, sigma_b^2), nrow = 2)
cafe_params <- mvrnorm(n_cafes, c(0, 0), Sigma)
a_cafe <- a_bar + cafe_params[, 1]
b_cafe <- b_bar + cafe_params[, 2]

# Simulate observations

n_per <- 10
d <- do.call(rbind, lapply(1:n_cafes, function(j) {
afternoon <- rep(c(0, 1), each = n_per / 2)
y <- a_cafe[j] + b_cafe[j] * afternoon + rnorm(n_per, 0, 0.5)
data.frame(cafe = factor(j), afternoon = afternoon, y = y)
}))

cat("Cafe data simulated:\n")
cat(" Cafes:", n_cafes, "\n")
cat(" Observations per cafe:", n_per, "\n")
cat(" True intercept SD:", sigma_a, "\n")
cat(" True slope SD:", sigma_b, "\n")
cat(" True correlation:", rho, "\n")

# Fit varying slopes model

library(lme4)
m <- lmer(y ~ afternoon + (1 + afternoon | cafe), data = d)

cat("\nVarying slopes model:\n")
cat(" Fixed: intercept =", round(fixef(m)[1], 3),
" afternoon =", round(fixef(m)[2], 3), "\n")
cat(" Random SD: intercept =", round(attr(VarCorr(m)$cafe, "stddev")[1], 3),
  "  slope =", round(attr(VarCorr(m)$cafe, "stddev")[2], 3), "\n")
cat(" Correlation:", round(attr(VarCorr(m)$cafe, "correlation")[1, 2], 3), "\n")

Click Run to execute this code in your browser with WebR

13.1.2 The LKJ prior for correlation matrices

The correlation matrix needs a prior. The LKJ distribution is a standard choice — it regularizes toward zero correlations:

$p(\mathbf{R}) \propto |\mathbf{R}|^{\eta - 1}$

When $\eta = 1$ , all correlation matrices are equally likely. When $\eta > 1$ , the prior concentrates toward $\rho = 0$ .

# Demonstrate LKJ prior for correlation matrices
# Using rethinking or manual simulation

set.seed(42)
n_samples <- 5000

# LKJ(eta = 1): uniform over correlations

# LKJ(eta = 2): favors rho near 0

# LKJ(eta = 4): strongly favors rho near 0

# Simple 2x2 case: sample correlations

# For 2x2, LKJ reduces to: p(rho) proportional to (1 - rho^2)^(eta - 1)

sample_lkj_2x2 <- function(eta, n) {

# Use Beta distribution trick

# rho ~ 2 * Beta(eta, eta) - 1

2 * rbeta(n, eta, eta) - 1
}

rho_1 <- sample_lkj_2x2(1, n_samples)
rho_2 <- sample_lkj_2x2(2, n_samples)
rho_4 <- sample_lkj_2x2(4, n_samples)

cat("LKJ prior on correlation (2x2 case):\n")
cat(" eta = 1: mean rho =", round(mean(rho_1), 3), " SD =", round(sd(rho_1), 3), "\n")
cat(" eta = 2: mean rho =", round(mean(rho_2), 3), " SD =", round(sd(rho_2), 3), "\n")
cat(" eta = 4: mean rho =", round(mean(rho_4), 3), " SD =", round(sd(rho_4), 3), "\n")

cat("\n Higher eta = more regularization toward rho = 0\n")

# Plot

par(mfrow = c(1, 3))
hist(rho_1, breaks = 50, col = "steelblue", main = "LKJ(eta=1)",
xlab = "rho", xlim = c(-1, 1))
hist(rho_2, breaks = 50, col = "steelblue", main = "LKJ(eta=2)",
xlab = "rho", xlim = c(-1, 1))
hist(rho_4, breaks = 50, col = "steelblue", main = "LKJ(eta=4)",
xlab = "rho", xlim = c(-1, 1))
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

Rethinking

Varying slopes are interactions. A varying slope is just an interaction between a predictor and a grouping factor — but with partial pooling. In a fixed-effects model, you would write y ~ x * group to get group-specific slopes. In a multilevel model, you write (1 + x | group) and the slopes are drawn from a common distribution. This means the slopes are shrunk toward the population average, with more shrinkage for groups with less data. The correlation between intercepts and slopes is the multilevel analog of allowing the interaction to covary with the main effect. It is all connected.

Summary of Section 13.1

Varying slopes allow the effect of a predictor to differ across groups.
The model estimates a covariance matrix with variances and a correlation $\rho$ .
Use (1 + x | group) in lme4 for varying intercepts and slopes.
The LKJ prior regularizes the correlation matrix toward $\rho = 0$ .
Varying slopes are pooled interactions — interactions with shrinkage.

Practice

Example: Admission decisions and gender