Longitudinal mediation analysis

Hello brms users,

In discussion with some colleagues today, we were wondering if you could estimate a time varying effect of a mediator using brms. I will have a play in the next couple of weeks but keen to know of anyone has tried or has ideas on wjat the model might look like.

Here is what I am thinking so far;
t = multiple time points
p = participants
m = mediator
z = outcome

# One implementation
z ~ (t+m|t/p),
t ~ m|t/p

#Another idea a colleague suggested
z ~ (t*m|p),
t ~ m|p

EDIT: I had more of a think and changed the nesting.

Does this model give the effect of the mediator and how it is changing over time ?

Oh disclaimer Im still new to this so I thought it was an interesting question, dont hold back if its dumb.

  • Operating System: Linux
  • brms Version: 2.3.1


If anyone wants to try the model Paul recommends below then I have uploaded some simulated data kindly provided by the authors in article linked below. ExampleYing.csv (9.2 KB)

To those that are interested, the conversation was spurred from this article where they implemented their own function.

Is the mediator also varying across time points?

1 Like

Yes in this case the mediator itself also varies across time (as well as its effect/coefficients).

Thanks heaps

I see. I am a little bit confused by you model formulas above. I see that your mediator is m, your outcome is z, but what is you initial exogenous variable? Or are we having different definition ins mind of what a mediator is?

Sorry it could very well be that they are wrong, they were just a punt in my mind as to how we could model the effect of m mediating t.

The dudes in the article linked above have written code to do dynamic mediation, where the mediator and its coefficients change across time. Im reasonably new to bayesian analysis but I was hoping there would be a way to implement this in brms and it would be way more flexible.

I have attached the DAG for the model they describe which probably explains better than I am right now.


The formula for this DAG is defined in the paper as;


Did this help clarify anything or do I need to hit the drawing board and think on this more ?

Good that we talked about that ;-)

In brms, you can try fitting the above described model as follows:

bf(y ~ s(t) + s(t, by = m) + s(t, by = x) + (1 | p)) +
  bf(m ~ s(t) + s(t, by = x) + (1 | p)) +

Wow, that is so awesome thanks so much !
I will try implementing this and report back in a few days.

EDIT: I added some simulated data that works with Paul’s formula above.

Would you in general suggest using splines for time-varying coefficients, or is there an alternativ model notation as well?

If we don’t specify the form of the time-varying coefficients and just say they are continuous functions, either splines or Gaussian processes are the way to go I think. At least in brms.

How would you write up the formula with Gaussian process?

Replace s() with gp().

Thanks! I’m not sure if I completely understood, but is either using s() or gp() the only way to model time-varying covariates? Or can you also use the “simple” coefficient in some way? And for a univariate model, would you also model time-varying covariates with s(..., by)?

it depends on how you want coefficients to vary over time but splines and GPs are perhaps the most flexible ones.

The approach works for univariate models as well.

Ok, one last question (hopefully), as I’m trying to better understand how to mix time-varying with non-time-varying covariates.

If I don’t want x to vary over time, I could write y ~ s(t) + s(t, by = m) + x? And another question, coming up after I have read the docs from ?mgcv::s. As smooths are centered, the main effect should also be included in the model:

# NOTE: smooths will be centered, so need to include fac in model....

Does this affect the above mentioned formula, so do I need to write something like

y ~ t + m + s(t) + s(t, by = m) + x

or does brms handle this differently?


If you read the doc of ?mgcv::s closely, you see that this centering only happens for categorical by variables, which we don’t have in the situation discussed here.

Why is a m = f(x) necessary? isn’t y ~ s(t) + s(t, by = m) + s(t, by = x) + (1 | p)) sufficient to get the unmediated X effect? Unless the effect of X on M is also required…

PS: the s(…, by ) structure works also with continuous variables?