I’m working with hierarchical models, but I’m having trouble using a different number of levels. I can’t use a tool like rstanarm, otherwise that would make things a lot easier (FWIW I’m using cmdstanpy)

Suppose I’m fitting the following:

```
outcomes ~ (1 | group1) + (1 | group2)
link = binomial
```

From what I understand this is roughly written in stan as:

```
// model_level_2.stan
data {
// ---- dimensions of things
int<lower=0> N;
int<lower=1> L_group1;
int<lower=1> L_group2;
// ---- data
int<lower=0> trials[N]; // initial trials
int<lower=0> outcomes[N]; // initial successes
int<lower=1, upper=L_group1> group1[N];
int<lower=1, upper=L_group2> group2[N];
}
parameters {
vector[L_group1] beta_group1;
vector[L_group2] beta_group2;
}
model {
// set your priors on beta_ ...
for (n in 1:N)
outcomes[i] ~ binomial_logit( trials[n], beta_group1[ group1[n] ] + beta_group2[ group2[n] ])
}
```

My question is: **What if I wanted to run this on an arbitrary number of groups?** For example, if I wanted to fit the following:

```
outcomes ~ (1 | group1) + (1 | group2) + (1 | group3)
link = binomial
```

Would I have to write *another* stan file called `model_level_3.stan`

where I’ve added terms for `group3`

and so on? That seems like a bad way to do it. Does stan support something like this, or do people have an approach for it?

Apologies if my models are wrong / bad. I’m open to feedback on that, but I feel like my question still stands 😇