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

Thanks!

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.

library(brms)

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)
summary(fit)

# you need the github version of brms for this to run
marginal_effects(fit)

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,
Jack

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).
step

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.
Thanks!

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
marginal_effects(fit3)

image

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'])

image

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

2 Likes

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.

Jack