另一种方法是将对 lmer 的调用包装在一个函数中，该函数将断点作为参数传递，然后根据断点使用优化最小化拟合模型的偏差。这最大化了断点的配置文件对数似然性，并且通常（即，不仅仅是这个问题）如果包装器内部的函数（在这种情况下为 lmer）找到最大似然估计值，条件是传递给它的参数，则整个过程找到所有参数的联合最大似然估计。

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)

要获得断点的置信区间，您可以使用配置文件可能性。例如，添加qchisq(0.95,1)到最小偏差（对于 95% 置信区间），然后搜索foo(x)等于计算值的点：

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

对于这个玩具问题，有些不对称，但精度还不错。如果您有足够的数据使引导程序可靠，则另一种方法是引导估计程序。

jbowman 提出的解决方案非常好，只是补充几点理论说明：

• 鉴于所使用的指标函数的不连续性，轮廓似然可能非常不稳定，具有多个局部最小值，因此通常的优化器可能不起作用。这种“阈值模型”的通常解决方案是使用更繁琐的网格搜索，评估每个可能实现的断点/阈值天的偏差（而不是介于两者之间的值，如代码中所做的那样）。请参阅底部的代码。

• 在这个估计断点的非标准模型中，偏差通常没有标准分布。通常使用更复杂的程序。参见下面的 Hansen (2000)。

• 引导程序在这方面并不总是一致的，请参见下面的 Yu（即将发表）。

• 最后，我不清楚你为什么要通过重新围绕 Days 来转换数据（即 bp - x 而不仅仅是 x）。我看到两个问题：

  1. 使用此程序，您可以创建人工天数，例如 6.1 天、4.1 天等。例如，我不确定如何解释 6.07 的结果，因为您只观察了第 6 天和第 7 天的值？（在标准断点模型中，6 到 7 之间的任何阈值值都应该给您相同的系数/偏差）
  2. b1 和 b2 的含义相反，因为 b1 天数减少，而 b2 天数增加？所以无断点的非正式测试是 b1 != - b2

这方面的标准参考是：

• 标准 OLS：Hansen (2000) 样本拆分和阈值估计，计量经济学，卷。68，第 3 期。（2000 年 5 月），第 575-603 页。
• 更奇特的模型：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 (即将出版) The Bootstrap in Threshold Regression”，计量经济学理论。

代码：

# 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)]

你可以试试MARS模型。但是，我不确定如何指定随机效果。 earth(Reaction~Days+Subject, sleepstudy)

这是一篇提出混合效应 MARS 的论文。正如@lockedoff 提到的，我在任何包中都没有看到任何相同的实现。