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.

@paul.buerkner, there’s a question for you at the bottom of this post :-)

I tried working more on this and made some progress. I wanted to extend @paul.buerkner s code by:

  • modeling intercept and slopes at a population-level so that only change points vary between participants.
  • obtaining more interpretable parameter estimates.

Simulated data

Here is some data with common slopes but varying-by-person change points:

# Data parameters
intercept = -2
slope1 = 1.0
slope2 = -0.5
breakpoints = runif(30, 2, 7)

# Predictors
df <- data.frame(
  age = rep(seq(0, 10, by=0.5), 30),
  person = factor(rep(1:30, each = 21))

# Response with per-person break point
df$y = intercept + ifelse(
    df$age <= changepoints[df$person],
    yes = slope1 * df$age,  # before change
    no = slope1 * changepoints[df$person] +
      slope2 * (df$age - changepoints[df$person])  # after change
df$y = rnorm(nrow(df), mean = df$y, sd = 0.5)

# Plot it
ggplot(df, aes(y=y, x=age, color=person)) + 
  geom_point() + 
  geom_line() + 
  theme(legend.position = 'none')

Specify model

We specify a slope-change only model (the intercept of slope2 at the the changepoint (change) is slope1 * change). The prior on change is important.

# The model
bform3 <- bf(
  y ~ Intercept + slope1 * age * step(change - age) +  # Section 1
    (slope1 * change + slope2 * (age - change)) * step(age - change),  # Section 2
  Intercept + slope1 + slope2 ~ 1,  # Fixed intercept and slopes
  change ~ 1 + (1|person),  # Per-person changepoints around pop mean
  nl = TRUE

# Priors
bprior3 <- prior(normal(0, 5), nlpar = "Intercept") +
  prior(normal(0, 2), nlpar = "slope1") +
  prior(normal(0, 2), nlpar = "slope2") +
  prior(uniform(0, 10), nlpar = "change")  # Within observed range

# Initial values
inits3 = list(list(
  slope1 = slope1,
  slope2 = slope2,
  Intercept = intercept

# Fit it!
fit3 <- brm(bform3, data = df, prior = bprior3, chains = 1, inits = inits3)

Inspect fit

The parameter recovery looks (very!) promising:

Group-Level Effects: 
~person (Number of levels: 30) 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(change_Intercept)     1.46      0.19     1.16     1.90 1.00      111      236

Population-Level Effects: 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept_Intercept    -2.02      0.05    -2.13    -1.92 1.01      778     1006
slope1_Intercept        1.02      0.02     0.98     1.06 1.01      828     1061
slope2_Intercept       -0.47      0.01    -0.50    -0.44 1.00     1079     1434
change_Intercept        4.31      0.27     3.74     4.85 1.01      160      259


The individual change points are in the correct magnitude (< |(7-2)/2| from the mean change point), but they seem not to capture the values used for simulation:

plot(breakpoints, ranef(fit3)$person[,1,'change_Intercept'])


Has the order of person been changed in ranef(fit3) or is there another simple explanation of this, @paul.buerkner?


I hope it is just the order. You can see the ranef order by looking at its column names. Since you defined person as a factor, the order might be 1, 10, 11, …, 19, 2, 20, etc. which might explain the problems you see.

This is not a bug but related to how R handles factors, but not sure if it is desriable behavior either.

Hi Paul,

Is using inv_logit here an identification requirement? I ask because I’d like to use a similar model, though my running variable runs from -100 to +100 and, thus, it doesn’t make much sense to constrain omega to be positive.

If it’s a constraint, I can rescale the data. I was just wondering why you’d used it here.


You could constrain omega to any range of values by scaling the inv_logit appropriately. In the example above I scaled it from 0 to 10 but you can also scale it from -100 to 100.

I’m also interested in fitting a linear mixed model with a random change point. However, my dataset has an additional categorical predictor, “group”, with 3 levels. I’m trying to incorporate a fixed effect of group into @paul.buerkner’s example code above, but am currently at a bit of a loss.

For reference, if we ignore the change point for a moment, my linear mixed model formula would look as follows:
y ~ age * group + (1 + age | person)

However, I have reason to believe that my data would be better fit by introducing a change point. Ideally, I’d like to estimate the effect of my categorical predictor “group” on the change point as well. As a start, I first tried adapting the code above by just including a main effect of “group” as follows:

bform <- bf(
    # intercept and main effect of group
    y ~ b0 + b1 * group +
        # pre-change slope
        b2 * (age - omega) * step(omega - age) +
        # post-change slope
        b3 * (age - omega) * step(age - omega),
    # intercept, slopes and change point varying by person
    b0 + b2 + b3 + omega ~ 1 + (1 | person),
    # fixed effect of group
    b1 ~ 1,
    nl = TRUE

However, this gave me the following error message:

Error: Factors with more than two levels are not allowed as covariates.

I saw that this issue came up previously, but I couldn’t really follow the solution in that thread. It seems I should incorporate the non-linear parameters in a linear formula, but I’m not sure how to go about this.

Thanks in advance for any advice!

I think you can also treat the beta parameters for the slope as linear models and add your group as predictors.

Then I think you’d specify the model as follows:

bform <- 
    y ~ b0 + b1 * (age - omega) * step(omega - age) +
    b2 * (age - omega) * step(age - omega),
    b0 + b1 + b2 + omega ~ 1 + group + (1 | person),
    nl = TRUE

Keep in mind that this is just off the top of my head and I haven’t fit a similar model myself. Others might know whether or not it makes sense.

P.S. The error you got just means that you need to recode your factor variable into a set of dummy variables.

1 Like