When you have poor diagnostics with complex models, it’s always recommend to start with the simplest model possible (e.g., only outcome) and slowly build up the model to see where the issues are introduced.
As I’ve mentioned in your previous threads, hierarchical parameters like:
target += normal_lpdf(theta4[j] | mu4,tau4);
Can result in divergences and so the non-centered parameterisation is a recommended alternative