I’m having a funny issue with a hierarchical model I’m fitting in Rstan. It runs fine, and the posterior predictive checks all seem to check out. The model is a negative binomial regression with three crossed random effects. I’ve used noncentered parameterization to avoid the funnel issue. However, the model only runs successfully without divergent transitions if I specify redundant mean values for each of the random intercepts.

I think that the more “correct” approach is supposed to be something like this:

```
model {
for(a in 1:AA){
a1[a] = b0 + sigma_a * raw_a[a];
}
for(b in 1:BB){
b2[b] = sigma_b * raw_b[b];
}
for(c in 1:CC){
c2[c] = sigma_c * raw_c[c];
}
for(i in 1:N){
log_y_hat = a1[AA[i]] + b2[BB[i]] + c2[CC[i]] + coef * x[i];
}
y ~ neg_binomial_2_log(log_y_hat, phi);
b0 ~ normal(0,1);
raw_a ~ normal(0,1);
raw_b ~ normal(0,1);
raw_c ~ normal(0,1);
sigma_a~exponential(1);
sigma_b~exponential(1);
sigma_c~exponential(1);
phi ~ gamma(1,2);
coef ~ lognormal(-1, 1);
}
```

But the model will only fit without divergent transitions if I include a redundant mean value for each of the random intercepts (a_mu, b_mu, c_mu)

```
model {
for(a in 1:AA){
a1[a] = a_mu + sigma_a * raw_a[a];
}
for(b in 1:BB){
b2[b] = b_mu + sigma_b * raw_b[b];
}
for(c in 1:CC){
c2[c] = c_mu + sigma_c * raw_c[c];
}
for(i in 1:N){
log_y_hat = b0 + a1[AA[i]] + b2[BB[i]] + c2[CC[i]] + coef * x[i];
}
y ~ neg_binomial_2_log(log_y_hat, phi);
b0 ~ normal(0,1);
raw_a ~ normal(0,1);
raw_b ~ normal(0,1);
raw_c ~ normal(0,1);
a_mu ~ normal(0,1);
b_mu~normal(0,1);
c_mu~normal(0,1);
sigma_a~exponential(1);
sigma_b~exponential(1);
sigma_c~exponential(1);
phi ~ gamma(1,2);
coef ~ lognormal(-1, 1);
}
```

Have any of you run into this before? Or is there something I’m missing? I’d appreciate your thoughts!