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

Question

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

Answer 1

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.

Answer 2

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[1] &
Days<search.range[2], "Days", drop=TRUE)))

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

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

Answer 3

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