Hi,

I have a model with a few levels to the data.

There are groups, students in groups, and tutors.

We are trying to model student scores.

There are approximately:

- 2500 groups
- 25000 students
- 1500 tutors

It is clear, from analyzing the data that there are very specific group means and group variances. Additionally students and mentors have a wide variety of ability.

I’ve built a basic model in Stan for this, with three intercepts (group, student, tutor) and a group-specific variance.

However it can take many hours for just 100 iterations, and the diagnostic reports that every step hit the maximum tree depth. I’ve tried increasing the tree depth to 50, but same results (I know 50 may be extreme, and was testing to see if that was the issue.)

Attached is the Stan file.

Once the intercepts and variances are working OK, I plan to add some other covariates (non-intercepts)

I’d love any suggestions that the group might have on how to improve this model.

```
data{
// Counts of things
int N;
int N_student;
int N_group;
int N_tutor;
// Dependent variable
real score[N];
// Intercepts
int<lower=1, upper=N_student> student[N];
int<lower=1, upper=N_tutor> tutor[N];
int<lower=1, upper=N_group> group[N];
}
parameters {
real mu_eta; // grand location
real mu_scale;
vector<lower=0>[N_group] g_sigma; // group specific scale
vector[N_group] g_eta; // group specific intercept
real g_scale; // scaling of group specific intercept
vector[N_student] s_eta; // student specific intercept
real s_scale; // scaling of student specific
vector[N_tutor] t_eta; // tutor specific intercept
real t_scale; // scaling of tutor specific
}
transformed parameters {
vector[N_group] g_re;
vector[N_student] s_re;
vector[N_tutor] t_re;
real mu;
mu = mu_scale * mu_eta;
g_re = g_scale * g_eta; // intercpet for group
s_re = s_scale * s_eta; // intercept for student
t_re = t_scale * t_eta; // intercept for tutor
}
model{
vector[N] theta;
vector[N] sigma;
// Priors
g_sigma ~ cauchy(0,2.5);
mu_eta ~ normal(0,1);
mu_scale ~ normal(0,1);
g_eta ~ normal(0,1);
g_scale ~ normal(0,1);
s_eta ~ normal(0,1);
s_scale ~ normal(0,1);
t_eta ~ normal(0,1);
t_scale ~ normal(0,1);
// model likelihood
theta = mu + g_re[group] + s_re[student] + t_re[tutor];
sigma = g_sigma[group];
score ~ normal(theta, sigma);
}
```