在 R 中拟合两个断点的分段回归：一条水平线，一条直线，然后是一条水平线 [关闭]答案

【问题标题】：piecewise regression with two breakpoints to be fitted in R: a horizontal line, a straight line and then again a horizontal line [closed]在 R 中拟合两个断点的分段回归：一条水平线，一条直线，然后是一条水平线 [关闭]
【发布时间】：2019-03-09 10:25:44
【问题描述】：

我想用 R 中的两个断点进行分段回归：首先是斜率为 0 的水平线，然后是直线，然后是斜率为 0 的水平线。这两个断点也应该拟合。

我的数据是这样的（我总共有 60 个类似的数据集）：

x <- c(1.306, 1.566, 1.736, 1.854, 2.082, 2.328, 2.650, 2.886, 3.162, 3.392) 
y <- c(176.4, 188.0, 193.8, 179.4, 134.4, 119.0, 66.2, 58.2, 58.2, 41.2)

有人知道怎么做吗？

【问题讨论】：

如果它可能有一些用处，我发现一个带有偏移量“a * exp(-0.5 * pow((x-b) / c, 2.0)) + Offset”的高斯峰值方程，带有参数a = 1.4554605225541448E+02, b = 1.5839534147665826E+00, c = -5.6899030064297440E-01, 偏移量 = 4.6871995528236809E+01 给出了 OK 拟合，R 平方 = 0.9865 和 RMSE

标签： r regression breakpoints piecewise

【解决方案1】：

使用nls 在x 上拟合一条具有最小值和最大值的线，如下所示。 a和b是交点的x值，.lin1是中间部分的截距，.lin2是中间部分的斜率。

fm <- nls(y ~ cbind(1, pmax(pmin(x, b), a)), alg = "plinear", start = list(a = 2, b = 3))

给予：

Nonlinear regression model
  model: y ~ cbind(1, pmax(pmin(x, b), a))
   data: parent.frame()
       a        b    .lin1    .lin2 
   1.774    2.764  425.463 -134.940 
 residual sum-of-squares: 530.7

Number of iterations to convergence: 5 
Achieved convergence tolerance: 6.489e-09

水平部分位于与交叉点的x值对应的y值处：

predict(fm, list(x = coef(fm)[1:2]))
## [1] 186.06667  52.53333

或者可以计算为对应于最小和最大x值的y值：

predict(fm, list(x = range(x)))
## [1] 186.06667  52.53333

我们可以像这样绘制点和拟合：

plot(y ~ x)
xx <- seq(min(x), max(x), length = 100)
p <- predict(fm, list(x = xx))
lines(p ~ xx, col = "red")

【讨论】：

这真是帮了大忙，谢谢！唯一的问题是我有 60 个这样的数据集，并且 x 的范围不同。因此，对于某些初始值有效，而对于某些无效。
例如另一个数据集是：x
我尝试了这些作为起始值：a0 = min(x, na.rm = FALSE) + 0.3; b0 = max(x, na.rm = FALSE) - 0.3
假设您已删除 NA，另一种可能性是 start = list(a = mean(head(x, k)), b = mean(tail(x, k)))，其中适当地选择了 k。您的起始值对应于 k=1。如果您的某些数据并未真正遵循此模型，那么您可能无法获得这些数据集的收敛性，因此期望完全自动的结果可能不现实。

【解决方案2】：

更新的答案需要更少的人际互动，仍然，基本上和以前的答案一样。

您可以将点粗略地分组为上和下自动点。一些过渡点可能与上点或下点分组，所以我们使用 boxplot.stats 消除任何看起来像这些组中的异常值。然后我们可以取平均值高低点来估计水平高度线。我们还使用非异常值的上下点来确定转换的 x 值。

HighLine = (2*max(y) + min(y))/3
HighPoints = which(y >= boxplot.stats(y[y>HighLine])$stats[1])
HighY = mean(y[HighPoints])

LowLine = (max(y) + 2*min(y))/3
LowPoints = which(y <= boxplot.stats(y[y<LowLine])$stats[5])
LowY = mean(y[LowPoints])

x1 = max(x[HighPoints])
x2 = min(x[LowPoints])

plot(x,y)
lines(c(min(x), x1,x2, max(x)), c(HighY, HighY, LowY, LowY))

【讨论】：