Now, I also suspect that there is still some unexplained variability in my data that can be different across the J groups.
Therefore, I would also like to insert in the model above for each group a different \sigma_j and then a hyper prior \tau_{\sigma}.
Iâ€™m not sure what the best way is to go about this and which distributions on \sigma are recommended.
I hope somebody here can me some guidance on this problem or pointers to relevant literature?

One approach that is used in psychology, searchable with the term â€ślocation-scale modelsâ€ť, is to a model the log variance in the same way as the mean is being modeled. I see here for code.

There are also some packages out there that can fit these for you.

LMMELSM: (Disclosure: My package) - Handles latent, multivariate mixed effects location scale models. It can also do univariate, observed MELSMs of course. Benefits: Can also model random effect variances, not just residual variance. Cons: Only handles one random grouping effect (e.g., no nested or crossed random effects currently; maybe when I have some time and motivationâ€¦)

Brms: Obviously, an amazing package. It can handle MELSMs and LSMs. Benefits: Handles nested/crossed random effects. Cons: Does not handle latent variances, or random effect variances.

Some others: ICCier (Also my package, but Iâ€™d recommend LMMELSM/brms over this). Nlme can do /some/ variance models, but the support isnâ€™t great. Hedeker has a couple of non-Bayesian MELSM packages too.