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.
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.
Best Answer
I may have found a solution to my problem, but could someone please confirm what I've done?
I slightly modified the above functions b1 and b2. I have called them b3 and b4:
The mixed model is similar, but b3 and b4 are used: