Estimating two break points in a broken stick model with random effects in R

Question

Five months ago, jbowman posted a very useful answer to estimate the break point in a broken stick model with random effects in R. I never use "computing" like ifelse and I would like to estimate two break points. I should write two other functions like b1 and b2 but I don't know how. Can someone please tell me how to do that in R? Thanks!

jbowman's code:

library(lme4)
str(sleepstudy)

#Basis functions
bp = 4
b1 <- function(x, bp) ifelse(x < bp, bp - x, 0)
b2 <- function(x, bp) ifelse(x < bp, 0, x - bp)

#Wrapper for Mixed effects model with variable break point
foo <- function(bp)
{
  mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy)
  deviance(mod)
}

search.range <- c(min(sleepstudy$Days)+0.5,max(sleepstudy$Days)-0.5)
foo.opt <- optimize(foo, interval = search.range)
bp <- foo.opt$minimum
bp
[1] 6.071932
mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy)

Answer 1

I've been trying to work this out as well (though not with random effects), and I can only get a little further along. I'd be happy if @jbowman could take a look.

A model that allows for a double break (i.e. with breaks at $k_1$ and $k_2$ with $k_2 \gt k_1$) with slopes that are continuous is

$ Y = \beta_0 +\beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \epsilon $

where

$$ X_1 = \left\{ \begin{array}{ll} X & if\ X \leq k_1\\ k_1 & if\ X \gt k_1\\ \end{array} \right. $$

$$ X_2 = \left\{ \begin{array}{ll} 0 & if\ X \leq k_1\\ X-k_1 & if\ k_1 \leq X \leq k_2\\ k_2-k_1 & if\ X \gt k_2\\ \end{array} \right. $$

$$ X_3 = \left\{ \begin{array}{ll} 0 & if\ X \leq k_2\\ X & if\ X \gt k_2\\ \end{array} \right. $$

To convert this to R code, I think you can do this:

X1 <- function(x, k1) ifelse(x<=k1, x, k1)
X2 <- function(x, k1, k2) ifelse(x<=k1, 0, ifelse(x<=k2, x-k1, k2-k1))
X3 <- function(x, k2) ifelse(x<=k2, 0, x)

And then for breakpoints bp1 and bp2

out.lm <- lm(y ~ X1(x, bp1) + X2(x, bp1, bp2) + X3(x, bp2), data=yourdata)
Y <- out.lm$coef[1] + out.lm$coef[2]*X1(x, bp1) + out.lm$coef[3]*X2(x, bp1, bp2) + out.lm$coef[4]*X3(x, bp2)

Sadly, when i apply this to my own problem it produces rubbish :( I'm no R programmer, hopefully someone with a bit more knowledge can help?

Ideally it'd be great to have a generic function to handle $n$ breakpoints.

Answer 2

Section 9.3 of this R code calculates random effect estimates for a broken stick model with multiple fixed break points.