Cross-Lagged Multilevel Panel Model


I’m referring to the paper “How to compare cross-lagged associations in a multilevel autoregressive model.” from Schuurman et al. ( The authors show how they use Bayesian multivariate response multilevel models (using Bugs, code available here: 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?


From quickly looking at the paper, I have two suggestions.

  1. You should use x1 and y1 as predictors of x2 and y2 (not the other way round). Same for the 2nd time point.
  2. You should set rescor to TRUE.


True, this is what I always do in SEM/lavaan - no idea, why I switched it in my notation above. Would you say that the above model represents a cross-lagged multilevel panel model? I’m not sure in which way the paper’s authors consider the n=1 autoregressive model… But maybe that’s already covered by using subject as random intercept (and more complex nested strucures would be written as nested random intercepts then…).


I do think it represents a cross-lagged panel model. From my understanding, the autoregressive part is just that y3 is regressed on y2 which is in turn regressed on y1 (same for x). But I haven’t read the paper thoroughly enough to be sure they mean exactly that.


Ok, thanks! I consider this question as solved, but probably do not tag the thread als “solved”, so others may get in into this discussion as well.


Just running one of these models myself trying to replicate Lavaan findings in BRMS.

Typically in a CLPM you model covariance in outcomes and predictors. Setting set_rescor(TRUE) addresses the former but do we need to do anything to capture the latter ?


CLPMs assume equal spacing between waves, but that isn’t always the case and I suspect that it isn’t necessary in a multilevel setting. Is it possible to extend the syntax to account for unequal spacings between waves? Maybe by including time as a covariate?


If you have covariates, you could include:

bf(x1 ~ covariate + (1 |id1|subject))
bf(y1 ~ covariate + (1 |id1|subject))


f1a <- bf(x1 ~ covariate + (1 |id1|subject))
f1b <- bf(y1 ~ covariate + (1 |id1|subject))
f2a <- bf(x2 ~ x1 + y1 + covariate + (1 |id1|subject))
f2b <- bf(y2 ~ x1 + y1 + covariate + (1 |id1|subject))
f3a <- bf(x3 ~ x2 + y2 + covariate + (1 |id1|subject))
f3b <- bf(y3 ~ x2 + y2 + covariate + (1 |id1|subject))

model <- brm(f1a + f1b + f2a + f2b + f3a + f3b + set_rescor(TRUE), data = data)

If I’m not mistaken, this should take care of the correlations for the initial predictors. The rest of the predictors are outcomes too so should be accounted for.


I guess we could indeed “weight” the autoregressive predictors according their distance or let the model find out itself (using an interaction or something similar) how the autoregressive effects change of distance between time points.