Skip to content
Statistical RethinkingStatistical Rethinking

12.3 More than one type of cluster

Real data often contain more than one clustering factor. Crossed and nested random effects extend the basic multilevel model to handle multiple sources of variation.

12.3.1 Crossed vs nested effects

  • Nested: Groups are hierarchical. Students within schools, schools within districts.
  • Crossed: Groups are not hierarchical. Actors and blocks, where each actor appears in each block.
# Multilevel model with crossed random effects
# Example: chimpanzee prosociality with actor and block effects

set.seed(42)
n_actors <- 7
n_blocks <- 6
n_trials <- 50

# Simulate crossed data

actor <- sample(1:n_actors, n_trials, replace = TRUE)
block <- sample(1:n_blocks, n_trials, replace = TRUE)
x <- sample(c(0, 1), n_trials, replace = TRUE) # prosocial on left

# True effects

a_bar <- 0
sigma_actor <- 1.5
sigma_block <- 0.5
a_actor <- rnorm(n_actors, 0, sigma_actor)
a_block <- rnorm(n_blocks, 0, sigma_block)
beta <- 0.5 # prosocial effect

# Generate outcome

logit_p <- a_bar + a_actor[actor] + a_block[block] + beta * x
p <- 1 / (1 + exp(-logit_p))
y <- rbinom(n_trials, 1, p)

d <- data.frame(actor = factor(actor), block = factor(block), x = x, y = y)

# Fit crossed multilevel model

library(lme4)
m <- glmer(y ~ x + (1 | actor) + (1 | block),
family = binomial, data = d)

cat("Crossed multilevel model:\n")
cat(" Fixed effects:\n")
cat(" Intercept:", round(fixef(m)[1], 3), "\n")
cat(" x (prosocial):", round(fixef(m)[2], 3), "\n")
cat(" Random effects SD:\n")
cat(" Actor:", round(attr(VarCorr(m)$actor, "stddev"), 3), "\n")
cat("    Block:", round(attr(VarCorr(m)$block, "stddev"), 3), "\n")

cat("\n Actor effects (deviations):\n")
actor_eff <- ranef(m)$actor[, 1]
for (i in 1:n_actors) {
cat(sprintf(" Actor %d: %.3f\n", i, actor_eff[i]))
}
Click Run to execute this code in your browser with WebR

12.3.2 Nested effects

When groups are nested (e.g., students within schools within districts), each level adds its own random effect:

μi=α+γdistrict[i]+γschool[i]+γstudent[i]\mu_i = \alpha + \gamma_{\text{district}[i]} + \gamma_{\text{school}[i]} + \gamma_{\text{student}[i]}

# Nested random effects: students within schools within districts
set.seed(42)
n_districts <- 5
n_schools_per <- 3
n_students_per <- 20

# Generate nested structure

district <- rep(1:n*districts, each = n_schools_per * n*students_per)
school <- rep(1:(n_districts * n_schools_per), each = n_students_per)
x <- runif(length(district), 0, 10)

# True effects

a_bar <- 2
sigma_district <- 1.0
sigma_school <- 0.5
sigma_student <- 0.3

a_district <- rnorm(n_districts, 0, sigma_district)
a_school <- rnorm(n_districts * n_schools_per, 0, sigma_school)

# Linear predictor

mu <- a_bar + a_district[district] + a_school[school] + 0.3 * x
y <- mu + rnorm(length(mu), 0, 1)

d <- data.frame(
district = factor(district),
school = factor(school),
x = x, y = y
)

# Fit nested model

library(lme4)
m <- lmer(y ~ x + (1 | district) + (1 | district:school), data = d)

cat("Nested multilevel model:\n")
cat(" Fixed effects:\n")
cat(" Intercept:", round(fixef(m)[1], 3), "\n")
cat(" x:", round(fixef(m)[2], 3), "\n")
cat(" Random effects SD:\n")
cat(" District:", round(attr(VarCorr(m)$district, "stddev"), 3), "\n")
cat(" School: ", round(attr(VarCorr(m)[[2]], "stddev"), 3), "\n")
cat(" Residual:", round(attr(VarCorr(m), "sc"), 3), "\n")
Click Run to execute this code in your browser with WebR
Rethinking

The structure of the data matters. Crossed and nested effects are not just statistical tricks — they reflect the structure of the data-generating process. If actors and blocks are crossed, you need a crossed model. If students are nested in schools, you need a nested model. Getting the structure wrong leads to biased estimates and incorrect uncertainty. The multilevel model forces you to think about how the data were collected and what grouping factors matter. This is not a burden — it is an opportunity to build a model that reflects the real world.

Summary of Section 12.3

  • Crossed effects: Groups are not hierarchical (e.g., actors × blocks).
  • Nested effects: Groups are hierarchical (e.g., students within schools).
  • Add random effects for each clustering factor: (1 | factor).
  • Use district:school notation for nested effects in lme4.
  • The model structure should reflect the data-generating process.