'brms' for piecewise glmm with random change point

Interfaces brms
1

I have a question about how I can use ‘brms’ for piecewise glmm with a random change point.
I would like to fit a piecewise glmm to model after how many cases a surgeon passed his/her learning curve (i.e. a plateau is reached). A simple model consists of 2 pieces of which the 2nd piece has no slope ( = 0) (see example figure attachedexample ). Since I have data from multiple surgeons I would also like the change point to be random.

I read the post on Piecewise Linear Mixed Models With a Random Change Point, however, I couldn’t figure out how to make the slope of the second piece zero.

Could you help me out? Many thanks in advance!

  • Operating System: macOS 10.14.5
  • brms Version: 2.9.0
2

Can you post your current work in progress model so that we can try to amend according to your needs?

3

Thanks Paul, I played a bit with some example data and got more or less what I had in mind.
The code I used:

 library(segmented)
 data(plant)

  library(brms)
  bform <- bf(
          y ~ b0 + b1 * (time - omega) * step(omega - time),
          b0 + b1 + alpha ~ 1 + (1|group),
          # to keep omega within the time range of 0 to 700
          nlf(omega ~ inv_logit(alpha) * 700),
          nl = TRUE
          )

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

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

  fit <- brm(bform, data = plant, prior = bprior, iter = 2000, chains = 1)
  summary(fit)

Plotting the results give me the following:
example_plant

Could you tell me why the change point appears to be gradual and not one point?
How do I get the change point for each ‘group’ from the model fit?

Many thanks in advance!
Gerjon

4

The change point appearing gradually may just be an artifact, but it may also be the result of the posterior uncertainty in the change point. Please also consider plotting the uncertainty in the plots.
You can extract the samples of alpha by using coef(fit, summary = FALSE)$alpha and then do the inv_logit transformation yourself to obtain omega (which is the changepoint, right?)

5

Thanks, that helps!