# Estimating the break point in a broken stick / piecewise linear model with random effects in R [code and output included]

Can someone please tell me how to have R estimate the break point in a piecewise linear model (as a fixed or random parameter), when I also need to estimate other random effects?

I've included a toy example below that fits a hockey stick / broken stick regression with random slope variances and a random y-intercept variance for a break point of 4. I want to estimate the break point instead of specifying it. It could be a random effect (preferable) or a fixed effect.

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)

#Mixed effects model with break point = 4
(mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy))

#Plot with break point = 4
xyplot(
Reaction ~ Days | Subject, sleepstudy, aspect = "xy",
layout = c(6,3), type = c("g", "p", "r"),
xlab = "Days of sleep deprivation",
ylab = "Average reaction time (ms)",
panel = function(x,y) {
panel.points(x,y)
panel.lmline(x,y)
pred <- predict(lm(y ~ b1(x, bp) + b2(x, bp)), newdata = data.frame(x = 0:9))
panel.lines(0:9, pred, lwd=1, lty=2, col="red")
}
)


Output:

Linear mixed model fit by REML
Formula: Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject)
Data: sleepstudy
AIC  BIC logLik deviance REMLdev
1751 1783 -865.6     1744    1731
Random effects:
Groups   Name         Variance Std.Dev. Corr
Subject  (Intercept)  1709.489 41.3460
b1(Days, bp)   90.238  9.4994  -0.797
b2(Days, bp)   59.348  7.7038   0.118 -0.008
Residual               563.030 23.7283
Number of obs: 180, groups: Subject, 18

Fixed effects:
Estimate Std. Error t value
(Intercept)   289.725     10.350  27.994
b1(Days, bp)   -8.781      2.721  -3.227
b2(Days, bp)   11.710      2.184   5.362

Correlation of Fixed Effects:
(Intr) b1(D,b
b1(Days,bp) -0.761
b2(Days,bp) -0.054  0.181


-

Another approach would be to wrap the call to lmer in a function that is passed the breakpoint as a parameter, then minimize the deviance of the fitted model conditional upon the breakpoint using optimize. This maximizes the profile log likelihood for the breakpoint, and, in general (i.e., not just for this problem) if the function interior to the wrapper (lmer in this case) finds maximum likelihood estimates conditional upon the parameter passed to it, the whole procedure finds the joint maximum likelihood estimates for all the parameters.

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)  To get a confidence interval for the breakpoint, you could use the profile likelihood. Add, e.g., qchisq(0.95,1) to the minimum deviance (for a 95% confidence interval) then search for points where foo(x) is equal to the calculated value: foo.root <- function(bp, tgt) { foo(bp) - tgt } tgt <- foo.opt$objective + qchisq(0.95,1)
lb95 <- uniroot(foo.root, lower=search.range[1], upper=bp, tgt=tgt)
ub95 <- uniroot(foo.root, lower=bp, upper=search.range[2], tgt=tgt)
lb95$root [1] 5.754051 ub95$root
[1] 6.923529


Somewhat asymmetric, but not bad precision for this toy problem. An alternative would be to bootstrap the estimation procedure, if you have enough data to make the bootstrap reliable.

-
Thank you -- that was very helpful. Is this technique called a two-stage estimation procedure, or does it have a standard name that I could refer to / look up? –  lockedoff Dec 13 '11 at 18:39
It's maximum likelihood, or would be if lmer maximized the likelihood (I think the default is actually REML, you need to pass a parameter REML=FALSE to lmer to get ML estimates). just estimated in a nested manner rather than all at once. I've added some clarification at the front of the answer. –  jbowman Dec 13 '11 at 19:40
I had some optimization problems and wide CIs when inverting the profile likelihood with my real data, but got narrower bootstrap CIs in my implementation. Were you envisioning a nonparametric bootstrap with sampling with replacement on subjects' data vectors? I.e., for the sleepstudy data, this would entail sampling with replacement from the 18 (subject) vectors of 10 data points, without doing any resampling within a subject's data vector. –  lockedoff Dec 14 '11 at 21:10
Yes, I was envisioning a nonparametric bootstrap as you describe, but partially that's because I don't know much about advanced bootstrap techniques that may (or may not) be applicable. The profile likelihood-based CIs and bootstrap are both asymptotically accurate, but it could well be that the bootstrap is significantly better for your sample. –  jbowman Dec 14 '11 at 23:56
+1 for the code –  Andrew Dec 15 '11 at 12:21
You could try a MARS model. However, I'm not sure how to specify random effects. earth(Reaction~Days+Subject, sleepstudy)