Example hierarchical models with 3+ levels of hierarchy?

Modeling
7

Hm, on actually inspecting the generated Stan code, I actually don’t think this gets at the kind of nested hierarchy with heterogeneity I was seeking to implement. The specification:

brms::bf(
          y ~ 1 + (1 | school) + (1 | district)
   ,  sigma ~ 1 + (1 | school) + (1 | district)
)

Allows for school and district random effects on the magnitude of variability among students, but still expresses that there is a single SD that conveys how schools vary from one another in their mean in each district. The idea I was trying to implement would have a separate SD per district to capture variability across districts in how schools’ means vary. Ditto separate SD per district to capture variability across districts in how schools’ SDs vary.

To be concrete, here’s some pseudocode of the generative process I’m thinking of (easiest to parse if you start at the innermost for loop):

for(this_district in districts){

	this_district_mean_school_mean = rnorm(1,mean_district_mean_school_mean,sd_district_mean_school_mean)
	this_district_sd_school_mean = exp(rnorm(1,mean_district_logsd_school_mean,sd_district_logsd_school_mean))

	this_district_mean_school_logsd = rnorm(1,mean_district_mean_school_logsd,sd_district_mean_school_logsd)
	this_district_sd_school_logsd = rnorm(1,mean_district_sd_school_logsd,sd_district_sd_school_logsd)

	for(this_school in schools_in_this_district){

		this_school_mean = rnorm(1,this_district_mean_school_mean,this_district_sd_school_mean)
		this_school_sd = exp(rnorm(1,this_district_mean_school_logsd,this_district_sd_school_logsd))

		for(this_student in students_in_this_school){
			score = rnorm(1,this_school_mean,this_school_sd)
		}
	}
}