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

Question

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

    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来拟合一种具有以下形式的矩形函数

任何人都可以帮助我或指出正确的方向吗？谢谢！

Answer 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)

Answer 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）以及其他参数（在对数尺度上，因为这是泊松模型的默认 link 函数） :

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

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

Broken stick (or piecewise) regression with 2 breakpoints

r

linear-regression

piecewise