Hello,

first of all, apologies if this has all been discussed elsewhere. If somebody could point me to a nice (accessible?) tutorial offering workable advice on how to deal with these issues, that would be much appreciated…

**The issue:**

I have a fairly complex (non-linear) multi-level model and I have a few Rhats > 1.01. Overall, the fit is pretty good, but I am assuming it would be nice to tighten these up a little before publishing.

From looking at the pairs plot + some intuition, the issue is that my fixed effects are generally negatively correlated with my random effects. For example, in the plot below, you can see some of the parameters in a pairs plot. We can see that the z parameter is negatively correlated with the corresponding fixed effect param.

This makes sense to me. In the model’s parameter blocks, I have:

```
parameters {
// These are all the parameters we want to fit to the data
////////////////////////////////////
// fixed effects
////////////////////////////////////
/* in order to allow for correlations between the
variables, these are all stored in a list
these include bA, bS (stick weight), and the two spatial
sigmas, along with the floor (chance of selectin an
item at random)
*/
array[K] real bA; // weights for class A compared to B
array[K] real b_stick; // stick-switch rates
array[K] real<lower = 0> rho_delta; // distance tuning
array[K] real rho_psi; // direction tuning
///////////////////////////////
// random effects
///////////////////////////////
array[K] vector[L] zA; // weights for class A compared to B
array[K] vector[L] z_stick; // stick-switch rates
array[K] vector[L] z_delta; // distance tuning
array[K] vector[L] z_psi; // direction tuning
real<lower = 0> sig_a;
real<lower = 0> sig_stick;
real<lower = 0> sig_delta;
real<lower = 0> sig_psi;
}
transformed parameters {
// combine fixed and random effects
array[K] vector[L] uA;
array[K] vector[L] u_stick;
array[K] vector[L] u_delta;
array[K] vector[L] u_psi;
for (kk in 1:K) {
uA[kk] = bA[kk] + sig_a*zA[kk];
u_stick[kk] = b_stick[kk] + sig_stick*z_stick[kk];
u_delta[kk] = rho_delta[kk] + sig_delta*z_delta[kk];
u_psi[kk] = rho_psi[kk] + sig_psi*z_psi[kk];
}
}
```

All of the `u_’ parameters have Rhat < 1.01. So great, the uncertainty in the b and u parameters cancels out to some extent?

I’m sure this must happen quite often with multi-level models.

**Questions:**

Are there any tricks to improve this?

At one point, I was using an LJK-prior and modeling the full covariance matrix. I removed this in the hope that a simplified model structure would decrease the time taken to fit the model!

Should I add it back in?

Thanks