# Advice/Peer Review: Piecewise Linear Model w/ Random Change Point and Discontinuity

**URL:** <https://discourse.mc-stan.org/t/advice-peer-review-piecewise-linear-model-w-random-change-point-and-discontinuity/11919>\
**Category:** brms\
**Created:** [November 14, 2019, 10:01pm UTC](https://discourse.mc-stan.org/t/advice-peer-review-piecewise-linear-model-w-random-change-point-and-discontinuity/11919 "2019-11-14T22:01:56Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jackbailey](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jackbailey/32/14978_2.png) [@jackbailey](https://discourse.mc-stan.org/u/jackbailey)\
**Post date:** [November 14, 2019, 10:01pm UTC](https://discourse.mc-stan.org/t/advice-peer-review-piecewise-linear-model-w-random-change-point-and-discontinuity/11919/1 "2019-11-14T22:01:56Z")

</div>

Operating System: macOS 10.14.6  
brms version: 2.10.0

Hi all,

[This popular post on piecewise linear models with random change points](https://discourse.mc-stan.org/t/piecewise-linear-mixed-models-with-a-random-change-point/5306/12) inspired me to try it out, but I’m having a little trouble fitting my model. In particular, I’d like to extend the model to include a discontinuity at the change point. Ultimately, I want to fit the following:

Y\_{i} \sim \beta\_{0} + \beta\_{1}(x\_{i} - \omega)I(x\_{i} \< \omega) + \beta\_{2}(x\_{i} - \omega)I(x\_{i} \geq \omega) + \beta\_{3}I(x\_{i} \geq \omega)

Where \beta\_{0} is the level of Y at the change point \omega, \beta\_{1} the slope before the change point, \beta\_{2} the slope after the change point, and \beta\_{3} the discontinuity at the change point. I(.) gives a value of 1 where its contents are true and 0 otherwise.

As you’ll see below, it doesn’t quite work yet. My suspicion is that the model isn’t identified. I can get the model to fit without the discontinuity, so that’s probably where the trouble is. But I don’t know why.

If anyone has any advice, I’d love to hear it.

**Example:** Below is code to reproduce this assuming a Gaussian outcome.

```
# Load packages
library(tidyverse)
library(brms)

# Create function to simulate data
sim_chgpt <- function(n = 100, b0 = 10, b1 = 0.5, b2 = 0, b3 = 30, om = 100){

# Create running variable for each observation
x <- runif(n, 0, 200)

# Add noise to omega parameter
om <- om + rnorm(n)

# Create outcome
y <-
  b0 +
  ( b1 * ( x - om ) * ( x < om ) ) +
  ( b2 * ( x - om ) * ( x >= om ) ) +
  ( b3 * ( x >= om ) ) +
  rnorm(n, 0, 5)

# Return tibble
tibble(x, y)
}

# Check graphically
sim_chgpt() %>% 
ggplot(aes(x = x, y = y)) +
geom_point()

# Fit model
m1 <-
  brm(
    formula = 
      bf(y ~
           b0 +
           b1 * ( x - omega ) * ( x < omega ) +
           b2 * ( x - omega ) * ( x >= omega ) +
           b3 * ( x >= omega ),
         b0 + b1 + b2 + b3 + omega ~ 1,
         nl = TRUE),
    family = gaussian(),
    prior = 
      prior(normal(0, 5), nlpar = "b0") +
      prior(normal(0, 5), nlpar = "b1") +
      prior(normal(0, 5), nlpar = "b2") +
      prior(normal(0, 20), nlpar = "b3") +
      prior(normal(100, 5), nlpar = "omega", class = "b"),
    data = sim_chgpt(n = 1000),
    iter = 2e3,
    refresh = 5,
    chains = 2,
    cores = 2
  )

```

This gives the following warnings:

```
Warning messages:
1: There were 1931 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
2: Examine the pairs() plot to diagnose sampling problems

3: The largest R-hat is 1.91, indicating chains have not mixed.
Running the chains for more iterations may help. See http://mc-stan.org/misc/warnings.html#r-hat 
4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See http://mc-stan.org/misc/warnings.html#bulk-ess 
5: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See http://mc-stan.org/misc/warnings.html#tail-ess 

```

And the following output:

```
Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b0_Intercept 2.00 8.20 -14.44 10.22 1.91 3 52
b1_Intercept -0.14 3.86 -10.01 7.88 1.85 95 70
b2_Intercept 0.24 0.23 -0.00 0.49 1.83 3 61
b3_Intercept 0.23 30.04 -38.11 30.87 1.84 3 94
omega_Intercept 49.99 49.88 -0.53 100.01 1.83 3 54

```

Clearly something is wrong. The estimates are off, the Rhats are high, and the bulk and tail ESSs are low.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [November 18, 2019, 2:19pm UTC](https://discourse.mc-stan.org/t/advice-peer-review-piecewise-linear-model-w-random-change-point-and-discontinuity/11919/2 "2019-11-18T14:19:08Z")

</div>

The discontinuity is IMHO indeed the problem: note that \beta\_1(x - \omega) I(x \< \omega) produces a likelihood continuous wrt. \omega, but \beta\_3(x \< \omega) does not. And Stan does not like non-continous likelihoods. To make this work, I would suggest you use `inv_logit` or similar to get a continous analog of the step function. This IMHO usually also makes more sense theoretically - I don’t think there are many processes where there is actual instantaneous change.

Hope that helps!

---

<div class="post-metadata">

**Author:** ![jackbailey](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jackbailey/32/14978_2.png) [@jackbailey](https://discourse.mc-stan.org/u/jackbailey)\
**Post date:** [November 18, 2019, 6:29pm UTC](https://discourse.mc-stan.org/t/advice-peer-review-piecewise-linear-model-w-random-change-point-and-discontinuity/11919/3 "2019-11-18T18:29:09Z")

</div>

Thanks for helping! It’s definitely the discontinuity: I fixed it to be positive and the problem went away (though there were a lot of max\_treedepth problems remaining).

Just so I’m clear, do you mean the following?:

```
formula = 
      bf(y ~
           b0 +
           b1 * ( x - omega ) * ( x < omega ) +
           b2 * ( x - omega ) * ( x >= omega ) +
           b3 * inv_logit( x >= omega ),
         b0 + b1 + b2 + b3 + omega ~ 1,
         nl = TRUE)

```

Or something else?

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [November 19, 2019, 5:49am UTC](https://discourse.mc-stan.org/t/advice-peer-review-piecewise-linear-model-w-random-change-point-and-discontinuity/11919/4 "2019-11-19T05:49:24Z")

</div>

I think `inv_logit( x >= omega )` won’t work as it is also discontinous, what I meant was `inv_logit( x - omega )`, which is continuous…

Also not sure how fixing `b3` to be positive would resolve your issues, so I consider it likely that the problem is sitll there (as witnessed by the max treedepth warnings).

Best of luck with your model!

---

<div class="post-metadata">

**Author:** ![jackbailey](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jackbailey/32/14978_2.png) [@jackbailey](https://discourse.mc-stan.org/u/jackbailey)\
**Post date:** [November 19, 2019, 5:02pm UTC](https://discourse.mc-stan.org/t/advice-peer-review-piecewise-linear-model-w-random-change-point-and-discontinuity/11919/5 "2019-11-19T17:02:20Z")

</div>

This seems to work, thanks!