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]))
}12.3.2 Nested effects
When groups are nested (e.g., students within schools within districts), each level adds its own random effect:
# 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")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:schoolnotation for nested effects inlme4. - The model structure should reflect the data-generating process.