@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)
```

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?