Piecewise regression for comparing drug switchers vs non-switchers

Hi!

Let’s start saying that I a newbie in the field of Bayesian Modelling, and just recently fell in love with the Bayesian framework thanks to @harrelfe after attending to one of his courses.
So apologies if my question is trivial

The setting: I want to compare some outcome y (e.g. lean muscle measured with dexa) in patients taking some medication A. Patients are followed up for years (not all with the exact same timing, but roughly once a year) and some of them might switch to drug B (actually the alternatives are more but let’s forget about it for now). So some patients will not switch during their fup, some will. And they switch not all at the same time (this is not a planned intervention at time t).

Let’s just thing for the moment at the simplest question: how can I compare the trends of y in switchers vs non-switchers? Most relevant the after switch slope compared to non-switchers?

I simulate some simple data with the following code.


# Parameters
intercept = -2
slope1 = 1.0
slope2 = -0.5
changepoints = runif(30, 2, 7)

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

df$y = rnorm(nrow(df), mean = df$z, sd = 0.5)
df$chpt = changepoints[df$person]

## Df for subjects without change

df_g <- data.frame(
  age = rep(seq(0, 10, by=0.5), 30),
  person = factor(rep(31:60, each = 21))
) 

df_g$z = intercept+1 + slope1 * df_g$age

## add some within subject noise
df_g$y = rnorm(nrow(df_g), mean = df_g$z, sd = 0.5)

## Merge them
df_g = df %>% mutate(Group="A") %>% 
  bind_rows(df_g %>% 
              mutate(Group="B",
                     chpt=0))

df_g = df_g %>% mutate(Group=factor(Group),
                       G=as.integer(Group=="A"))

## Plot
ggplot(df_g, aes(y=y, x=age,colour = Group)) + 
  geom_point() + 
  theme(legend.position = 'none') + 
  geom_smooth()

So I guess I would like to estimate

  • a common slope before switch. This slope is likely to be the only one for non-switchers
  • a slope after switch (for switchers only)
  • a population average time at switch.

I am assuming there is a substantial change after switch, but there might be little in reality

I read the post Piecewise Linear Mixed Models With a Random Change Point but here I know when the change point for each patient (among the switchers).

I modelled this in the past (traditionally) setting to zero the switch time for all switchers and median switch time for non switchers (also did with switcher-non switcher matching on time but not happy about it).

Any hint?

  • Operating System: Mac OS and Linux Debian
  • brms Version: 2.22.0

I don’t actually know the brms syntax, but I’m guessing that what you need to do is introduce a covariate for before/after change point. Then you can model a difference in slope before or after the switch. In Stan, that’d just be something like this:

parameters {
  real slope;
  real diff;
  ...
}
model {
  for (n in 1:N) {
    for (t in 1:T) {
      real adj_slope = t > change[n] ? slope + diff : slope;
      ... use adjusted slope ...
    }
  }

In brms, I’m guessing you can use an indicator covariate change[n, t] for subject n at time t that is 1 if t is after subject n’s changepoint. Then you can include the original slope and a new slope for the difference using covariate change.

2 Likes

Thanks @Bob_Carpenter

As I said I pretty new to Bayesian modelling and have a very vague idea of Stan syntax. I guess it would be worth investing some time at learning it.

Back to the point.

I actually did something alike with two slopes, slope_1 (before, if t < change) and slope_2 (after, i.e if t >= change) ).

First issue is that I have different groups, so ideally two slopes each.
Second, and more tricky, one group actually never switched, so there is no change-time.

What’s the trend for those who switch A->B? Those who switched A->C? And those who didn’t switch at all?

I currently cannot access the stan forums directly so can only respond via email. From what you describe, I think it would make most sense if you directly specified your model in Stan, building on the code that Bob provided. Your model may be possible in brms, but it would likely be more cumbersome to write in brms than in Stan directly.

I’m sure we can sort this out. Do you know why? Do you need us to generate another login or something for you? (If you can send email, I’ll see it—I may not see a reply to this forum post).

1 Like

It now seems to work again. It was some temporary issue with discourse perhaps.