I’m referring to the paper “How to compare cross-lagged associations in a multilevel autoregressive model.” from Schuurman et al. (http://dx.doi.org/10.1037/met0000062). The authors show how they use Bayesian multivariate response multilevel models (using Bugs, code available here: http://supp.apa.org/psycarticles/supplemental/met0000062/met0000062_supp.html) to fit a cross-lagged panel model in a multilevel design.

My question is, how this would translate into brms, and if this is currently possible? I have no such comprehensive statistical background to understand the formulas in detail. What I find confusing is that the authors state at p. 208 that *“The model as defined in Equations 1–3 forms level one of the multilevel model.”* - which seems to me that level 1 includes all observations. However, some sentences later they say that *“…whereas the model defined at level one is identical to an n=1 autoregressive time series model”*. So is level 1 defined as a n=1 model, or does it include all observations, or am I confusing something here?

Let’s say, we have 3 time points and **x** and **y** were measured at each time point, while we have **covariates** that are constant accross all time points. From what I understood, the model would be written in brms like this:

*Edited Code Example, this addresses Paul’s first comment*

```
f1a <- bf(x2 ~ x1 + y1 + covariate + (1 |id1|subject))
f1b <- bf(y2 ~ x1 + y1 + covariate + (1 |id1|subject))
f2a <- bf(x3 ~ x2 + y2 + covariate + (1 |id1|subject))
f2b <- bf(y3 ~ x2 + y2 + covariate + (1 |id1|subject))
# maybe also set_rescor(TRUE)?
model <- brm(f1a + f1b + f2a + f2b + set_rescor(FALSE), data = data)
```

Is this correct? Would this describe a cross-lagged panel model with random intercept structure?