我是在这里提问的新手，所以我希望这是一个合理的地方。
我正在尝试拟合分段回归。我希望我的响应 (y) 从 x0 增加到沿 x 的某个位置，然后在断点后达到稳定；此后，斜率=0。我想拟合的模型有 3 个参数——截距、斜率和断点。在最佳情况下，我会将此模型与使用 AIC 或类似方法的仅截距/空模型进行比较。
我从包分段开始，这个提示将第二个斜率修复为 0： https ://stat.ethz.ch/pipermail/r-help/2007-July/137625.html 但是在更新了一堆包之后，整个事情（已经有点挑剔了）完全停止了工作。
所以，现在，我正在使用 optim()，但我不确定如何正确估计错误，尤其是在断点参数附近。有关使用 bblme 和 nls 的尝试示例，请参见下面的代码（它们都有效，但在我的数据集上似乎都失败了）。
根据一位使用 SAS 的同事的说法，proc nlin 似乎能够轻松地拟合这个模型并为参数生成标准误差估计......但是，我决定不使用 SAS。
这种交流（和其他交流）似乎表明人们普遍担心这种模型的错误估计，特别是断点参数 - 它引用了 Venables 的一篇可能有用的帖子，不幸的是不再可用。 使用 R 的 nls() 进行变化点分析
（注意，我的数据由一个代表努力的变量加权。当它可用时，您可以将 wt 提供给 weights=() 语句，否则，我使用扩展的数据集或加权模型中的平方和。）
模拟样本量减少的分段数据
set.seed(15)
a.sim<-0 # intercept
b.sim<-0.5 # slope for segment 1
n<-12
brk<-4 # breakpoint
x <- c(1:n)
y <- numeric(n)
y[1:brk]<-a.sim+b.sim*x[1:brk]+rnorm(brk, 0, .2)
y[(brk+1):n] <- (a.sim+ b.sim* brk) + rnorm((n-brk), 0, .2) # second, flat segment
y[n]<-y[n]-.30*y[n] #subtract from last point, as these are often outside of the normal pattern
wt<-c(rep(50, n-4), c(40, 40, 35, 5)) #weight variable
dat<-as.data.frame(cbind(x, y, wt)) # make dataframe
dat.expand <- dat[ rep( seq(dim(dat)[1]), dat$wt),]# second dataset with rows repeated based on weight
with(dat, plot(x,y, ylim=c(0, max(y)), pch=16, cex=wt/(10)))### plot, with symbols representing weight variable
使用 optim 解决
mod<-function(par, data){
a <- par[1]
b <-par[2]
x.star <-par[3]
yfit<-c(NA,length(data$y))
small.x<-I(data$x<=x.star)
yfit[small.x==T]<-(a+b*data$x[small.x==T])
yfit[small.x!=T]<-(a+b*x.star)
sum((data$y-yfit)^2)
}
fit1<-optim(par=c(a=.5, b=.5, x.star=3), mod, data=dat, hessian=T)
这几乎总是能很好地工作并提供合理的配合。对于生产运行，我给它提供了一个潜在起始值的大表，以便在实际优化之前找到一个合理的起始点。
独联体的 bblme
我与一位统计学家交谈过，他建议将此作为在断点周围生成置信区间的一种可能方法。在某些情况下它适用于我的数据，但通常会出错。
library(bbmle)
mod2<-function(a,b,x.star,sig){
yfit<-c(NA,length(y))
small.x<-I(x<=x.star)
yfit[small.x==T]<-(a+b*x[small.x==T])
yfit[small.x!=T]<-(a+b*x.star)
-sum(dnorm(y,yfit,1, log=TRUE))}
fit3<-mle2(mod2, start=list(a=0,b=0.5,x.star=brk), data=dat)
ci<-confint(fit3)
方法
这基本上不适用于我的数据类型（如果我模拟更大的数据集则有效）我得到奇异梯度矩阵误差或......步长因子 0.000488281 降低到 0.000976562 的“minFactor”以下......
nls.mod<-nls(y ~ ifelse(x <x.star, a+b*x,a+b*x.star),
data = dat, weights=wt,
start = c(x.star=brk,b=0.4, a=0))
JAGS？
所以，在统计学家的帮助下，我也在 JAGS/r2jags 中进行了这项工作。如果有一种直接的、频繁的方法来合理估计断点周围的错误，我宁愿不使用这种方法。我不是在这里寻求完美。只是一些合理的东西..