Piecewise Linear Mixed Models With a Random Change Point

Interfaces brms
9

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

image
image.png700×432 23.7 KB

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?

5 Likes