Example hierarchical models with 3+ levels of hierarchy?

Modeling
2

I just coded up a quick example using brms. Does this help you?


# number of students/scores
N_g <- 50
# number of schools
N_j <- 20
# number of districts
N_k <- 10
# total cases
N <- N_g * N_j * N_k

# intercept (e.g. IQ scores)
alpha <- 100
student_vec <- 1:N
school_vec <- rep(1:N_j, each=N_g, times = N_k)
district_vec <- rep(1:N_k, each=N_j*N_g)
# random intercept for schools
r_j <- rnorm(N_j, 0, 10)
# random intercepts for districts
r_k <- rnorm(N_k, 0, 5)
# linear predictor
mu <- alpha + r_j[school_vec] + r_k[district_vec]
# noise
sigma <- 5

# simulate data
y <- rnorm(N, mu, sigma)

# create data frame
d <- data.frame(
  y=y, 
  school = school_vec, 
  district = district_vec
)

# fit the model (can take a few minutes depending on N)
fit <- brms::brm(
  y ~ 1 + (1 | school) + (1 | district), 
  data=d, 
  cores=4
)

summary(fit)

And the output from a representative run:

Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y ~ 1 + (1 | school) + (1 | district) 
   Data: d (Number of observations: 10000) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~district (Number of levels: 10) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     3.72      1.06     2.29     6.17        622 1.01

~school (Number of levels: 20) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)    11.44      1.91     8.25    15.86        512 1.00

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept   100.73      2.81    95.08   106.36        369 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     4.99      0.03     4.93     5.06       1389 1.00
2 Likes