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

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)  - Are you interested in having the slopes be continuous at the breakpoints or not? – jbowman May 29 '12 at 13:14 add comment ## 2 Answers 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.

-
A special note to get @jbowman 's attention as I can't add a comment to the question (too little rep). –  a different ben Aug 29 '12 at 7:05
Sorry, been awfully busy over the last few months, and I'll take a look at this over the next few days (whether or not you're interested any more!) –  jbowman Dec 27 '12 at 17:20
Yes indeed @jbowman, still interested! –  a different ben Jan 11 at 4:24

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

-