带有 2 个断点的断棒（或分段）回归答案

【问题标题】：Broken stick (or piecewise) regression with 2 breakpoints带有 2 个断点的断棒（或分段）回归
【发布时间】：2018-09-24 13:50:55
【问题描述】：

我想用下一个数据估计一个函数的两个断点：

    df = data.frame (x = 1:180,
                y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 2, 2, 4, 2, 2, 3, 2, 1, 2,0, 1, 0, 1, 4, 0, 1, 2, 3, 1, 1, 1, 0, 2, 0, 3,  2, 1, 1, 1, 1, 5, 4, 2, 1, 0, 2, 1, 1, 2, 0, 0, 2, 2, 1, 1, 1, 0, 0, 0, 0, 
                    2, 3, 0, 3, 2, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
# plotting y ~ x 
plot(df)

我知道该函数有两个断点：

y = y1 if x < b1;
y = y2 if b1 < x < b2;
y = y3 if b2 < x;

而我想找到b1和b2来拟合一种如下形式的矩形函数

谁能帮助我或指出正确的方向？谢谢！

【问题讨论】：

如果你能告诉我们数据来自哪里，那么很可能已经有一些特定领域的 R 包可以解决这个问题。
@zx8754 来源于动物行为实验心理学； y 是响应率 (responses/sec) 作为经过时间 (x, in sec) 的函数。我想使用简单的 OLS。

标签： r linear-regression piecewise

【解决方案1】：

1) kmeans 像这样尝试kmeans：

set.seed(123)
km <- kmeans(df, 3, nstart = 25)

> fitted(km, "classes") # or equivalently km$cluster
  [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [38] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

> unique(fitted(km, "centers")) # or except for order km$centers
      x         y
3  30.5 0.5166667
1  90.5 0.9000000
2 150.5 0.0000000

> # groups are x = 1-60, 61-120 and 121-180
> simplify2array(tapply(df$x, km$cluster, range))
       1   2  3
[1,]  61 121  1
[2,] 120 180 60

plot(df, col = km$cluster)
lines(fitted(km)[, "y"] ~ x, df)

2) 蛮力法 另一种方法是蛮力法，我们计算每一对可能的断点，并选择线性模型中平方和最小的对。

grid <- subset(expand.grid(b1 = 1:180, b2 = 1:80), b1 < b2)

# the groups are [1, b1], (b1, b2], (b2, Inf)
fit <- function(b1, b2, x, y) {
   grp <- factor((x > b1) + (x > b2))
   lm(y ~ grp)
}

dv <- function(...) deviance(fit(...))

wx <- which.min(mapply(dv, grid$b1, grid$b2, MoreArgs = df))

grid[wx, ]
##       b1 b2
## 14264 44 80

plot(df)
lines(fitted(fit(grid$b1[wx], grid$b2[wx], x, y)) ~ x, df)

【讨论】：

有趣的方法。两个问题：1）我怎样才能知道方法与数据集的拟合程度？ 2）这种方法对相同值的点数和点的值都敏感吗？
(1) 您可以检查平方和之间的平方和占总平方和的百分比。只需输入km 即可查看。 (2) 它试图最小化残差的平方和，其中一个点的残差是到其集群中心的距离。
添加了第二种（非等效）方法。
第二个对我来说似乎更自然。谢谢！
关于2）b1的范围比b2宽有什么原因吗？

【解决方案2】：

我可以看到 y 是整数，所以最好使用泊松或二项式模型来估计。这是使用R包mcp的解决方案：

# Three intercept segments
model = list(
  y ~ 1,
  ~ 1,
  ~ 1
)

library(mcp)
fit = mcp(model, df, family = poisson(), par_x = "x", adapt = 2000)
plot(fit)

请注意，mcp 是唯一可以估计变化点 ant 参数估计周围的不确定性的软件包之一。摘要显示了估计变化点的位置（cp_1 和 cp_2）以及其他参数（在对数尺度上，因为这是泊松模型的默认链接函数）：

summary(fit)

Family: poisson(link = 'log')
Iterations: 9000 from 3 chains.
Segments:
  1: y ~ 1
  2: y ~ 1 ~ 1
  3: y ~ 1 ~ 1

Population-level parameters:
  name   mean lower  upper Rhat n.eff
  cp_1  39.57  37.8  45.00    1    54
  cp_2  99.82  99.0 101.21    1  2211
 int_1  -4.00  -6.5  -1.88    1   577
 int_2   0.32   0.1   0.54    1  6288
 int_3 -11.02 -20.9  -3.56    1  2487

【讨论】：