# 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 • Any way to make bp a random effect? Feb 25 '14 at 21:42

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  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, upper=bp, tgt=tgt)
ub95 <- uniroot(foo.root, lower=bp, upper=search.range, tgt=tgt)
lb95$root  5.754051 ub95$root
 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? 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. 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. 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. Dec 14 '11 at 23:56

The solution proposed by jbowman is very good, just adding a few theoretical remarks:

• Given the discontinuity of the indicator function used, the profile-likelihood might be highly erratic, with multiple local minima, so usual optimizers might not work. The usual solution for such "threshold models" is to use instead the more cumbersome grid search, evaluating the deviance at each possible realized breakpoint/threshold days (and not at values in between, as done in the code). See code at bottom.

• Within this non-standard model, where the breakpoint is estimated, the deviance does usually not have the standard distribution. More complicated procedures are usually used. See the reference to Hansen (2000) below.

• The bootstrap is neither always consistent in this regard, see Yu (forthcoming) below.

• Finally, it is not clear to me why you are transforming the data by re-centering around the Days (i.e., bp - x instead of just x). I see two issues:

1. With this procedure, you create artificial days such as 6.1 days, 4.1 etc. I am not sure how to interpret the result of 6.07 for example, since you only observed values for day 6 and day 7 ? (in a standard breakpoint model, any value of the threshold between 6 and 7 should give you same coef/deviance)
2. b1 and b2 have the opposite meaning, since for b1 days are decreasing, while increasing for b2? So the informal test of no breakpoint is b1 != - b2

Standard references for this are:

• Standard OLS: Hansen (2000) Sample Splitting and Threshold Estimation, Econometrica, Vol. 68, No. 3. (May, 2000), pp. 575-603.
• More exotic models: Lee, Seo, Shin (2011) Testing for threshold effects in regression models, Journal of the American Statistical Association (Theory and Methods) (2011), 106, 220-231
• Ping Yu (forthcoming) The Bootstrap in Threshold Regression", Econometric Theory.

Code:

# Using grid search over existing values:
search.grid <- sort(unique(subset(sleepstudy, Days > search.range &
Days<search.range, "Days", drop=TRUE)))

res <- unlist(lapply(as.list(search.grid), foo))

plot(search.grid, res, type="l")
bp_grid <- search.grid[which.min(res)]


You could try a MARS model. However, I'm not sure how to specify random effects. earth(Reaction~Days+Subject, sleepstudy)

• Thanks -- I browsed through the package documentation but it did not seem to support random effects. Dec 13 '11 at 18:39

This is a paper that proposes a mixed effects MARS. As @lockedoff mentioned, I don't see any implementations of the same in any package.