Piecewise Linear Mixed Models With a Random Change Point

I’m wondering if it’s possible to implement a piecewise linear mixed model with a random change point in brms? This paper implements the model in Stan:

Bayesian Piecewise Linear Mixed Models With a Random Change Point

  • Operating System: Windows
  • brms Version: 2.4.0


I haven’t read the paper in detail, but you may try the following model. I built it rather ad-hoc so not sure if that’s the best solution we have in brms. Please also note that the priors are just to get the model running with this non-sensible simulated data.


bform <- bf(
  y ~ b0 + b1 * (age - omega) * step(omega - age) + 
    b2 * (age - omega) * step(age - omega),
  b0 + b1 + b2 + alpha ~ 1 + (1|person),
  # to keep omega within the age range of 0 to 10
  nlf(omega ~ inv_logit(alpha) * 10),
  nl = TRUE

df <- data.frame(
  y = rnorm(330),
  age = rep(0:10, 30),
  person = rep(1:30, each = 11)

bprior <- prior(normal(0, 3), nlpar = "b0") +
  prior(normal(0, 3), nlpar = "b1") +
  prior(normal(0, 3), nlpar = "b2") +
  prior(normal(0, 3), nlpar = "alpha")

make_stancode(bform, data = df, prior = bprior)

fit <- brm(bform, data = df, prior = bprior, chains = 1)

# you need the github version of brms for this to run

This looks cool. I’ll try it out. Thanks, Paul

Nice piece of code. I am intrigued about what the step function is doing in this instance- can anyone shed any light on it?

Similarly, how would one develop the code above to allow for a quadratic relationship in either/both slopes?

thanks in advance,

Stumbled upon this when I was thinking if there’s anything brms couldn’t do.

Actually I learned quite a bit of syntax that was completely new to me, e.g. add up response variables is just a shorthand way of giving them the same link functions.

for the question on step function, I tried the RStudio help (which direct to something related to AIC), and googling “r formula step”, “brms formula step”. No luck.

Then it dawned on me that brms relies on stan to do the heavy-lifting of HMC, so maybe it is a stan function? then I actually tried to read the code generated by the make_stancode function as suggested by @paul.buerkner at #2 post. And there it was!

So I go back to google “stan function step”, and here’s the result (you have to scroll down a bit).

And it does a simple thing, if the parameter is positive, it returns 1, otherwise 0.

1 Like

Dear Paul,
thank you for outlining the model formula above, which is really helpful! I have a rich dataset with a mean of 100 repeated observations per individual and reason to suspect that there is actually more than one random change point. Is it possible to add a third slope (b3) and a second random change point (between b2 and b3), and how would one go about it? Particularly, I am not sure how to add a second alpha term.

I am not sure exactly, but one way to do it would be to introduce a new omega (and new alpha if needed) and to add them to the model. Then you need to make sure somehow that one of the omegas is always larger then the other to identify the model.

In any case, it may make sense that you first write down the mathematical model you want to estimate so that we have a basis for thinking of how to write that down in brms or Stan.

Thanks for the quick answer. I see if I can come up with mathematical model.