Apparently equivalent noncentered parameterization, different outcomes

Modeling
,
1

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!

2

Why do you feel that these parameters are redundant? In one case, you have a gaussian multiplied by an exponential, and, in the other, a mixture of a gaussian and a gaussian mutliplied by an exponential. Their geometries are not the same:

3

Hi Corey,
Thanks, I see what you mean. I was thinking in terms of the intercepts’ role in the model rather than their geometries. Since one intercept is already specified (b0), I thought that the random intercepts would only need to be included as some deviation from that mean value. I don’t think it’s a problem to have those three extra intercept terms, but I’m hoping to understand why they are necessary for the model fit. (Right now I think my response to any reviewer questions about it would be kind of lame.)