I am creating a gentle introduction to seemingly unrelated regression (SUR) models for a class of novices. I was looking into rstanarm’s function for multivariate outcomes, stan_mvmer. However, it seems that that this function does not model the correlation between the residuals. Is that right? (Although I think it would be a useful exercise for students to calculate residuals and correlations themselves, that’s probably too much for a first class session.)

If not rstanarm, are there other recent scripts for multivariate outcomes? My inclination is to start with convenience packages like rstanarm, if possible. But the class will have some familiarity with stand-alone Stan models by this point, so working with Stan models directly is also an option, if necessary.

Here’s a related thread about what’s not currently implemented in stan_mvmer:

Like @tjmahr said, brms has a bunch of stuff you can do with multivariate outcomes. It doesn’t do everything that @sambrilleman did with stan_jm and stan_mvmer in rstanarm, but it also does some things that can’t be done with those rstanarm functions.

Actually I find the example Stan code listed in section 9.15. of the current manual version also not that complicated and well suited for educative purposes. So why not use it (you might skip the LKJ part for the sake of clarity)?