Snijders & Bosker (2012, chapter 16) suggest using the following syntax in `nlme`

for fitting a multivariate multilevel model (data & scripts for the chapter can be found here)

```
library(nlme)
mlb161 <- lme(response ~ - 1 + whichpart, random = ~ -1 + whichpart|schoolnr,
weights=varIdent(form=~1|whichpart),
corr=corSymm(form=~as.numeric(whichpart)|schoolnr/pupilNR_new),
data=mlbook_disaggregated, method="ML",
control = list(maxIter=500, msMaxIter=500, tolerance=1e-8,
niterEM=250))
```

This assumes that the residuas are clustered within `pupilNR_new`

(which represents students nested in schools, i.e., `schoolnr`

): `corr=corSymm(form=~as.numeric(whichpart)|schoolnr/pupilNR_new)`

Taking this model to a Bayesian framework and fitting it with `brms`

, I would like to know what would be the equivalent of this type of residuals clustering. I understand that the `autocor`

or `rescor`

arguments should be relevant for defining it, but I’m not sure how it could be specified in an homologous way as defined with the `nlme::lme()`

example using `corSymm()`

.

In sum, taking aside the change from a frequentist to a bayesian perspective, what would be the proper way to specify this multivariate multilevel model in `nlme`

with `brms`

? (I’ve also made this question in Stack Overflow with no luck)

**References**

Snijders, T. A. B., & Bosker, R. J. (2012). *Multilevel analysis: An introduction to basic and advanced multilevel modeling* (2nd ed). Sage.