分段回归，包含一条直线和一条水平线在分界点相交

Question

我希望进行一次有一个拐点的分段线性回归，其中回归线的第二部分的斜率为0。有关如何进行分段线性回归的示例，例如这里。我遇到的问题是我不清楚如何将模型的一半的斜率固定为0。

我尝试过：

lhs <- function(x) ifelse(x < k, k-x, 0)
rhs <- function(x) ifelse(x < k, 0, x-k)
fit <- lm(y ~ lhs(x) + rhs(x)) 

其中k是断点，但右侧的线段不是平坦/水平的。

我想将第二个线段的斜率限制为0。我尝试了：

fit <- lm(y ~ x * (x < k) + x * (x > k))

但是，我不确定如何使第二半段具有零斜率。

非常感谢任何帮助。

---

我的解决方案

由于下面的评论，我找到了一个解决方案。以下是我用来优化和绘制拟合曲线的代码:

x <- c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1)
y <- c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516, 
0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527, 
0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.25414608)

range_x <- max(x) - min(x)
intervals <- 1000
coef1 <- c()
coef2 <- c()
r2 <- c()

for (i in 1:intervals) {
  k <- min(x) + (i-1) * (range_x / intervals)     
  x2 = (x - k) * (x < k)
  fit <- lm(y ~ x2)
  coef1[i] <- summary(fit)$coef[1]
  coef2[i] <- summary(fit)$coef[2]
  r2[i] <- summary(fit)$r.squared
  }

best_r2 <- max(r2)   # get best r squared
pos <- which.max(r2)                                          
best_k <- min(x) + (pos - 1) * (range_x / intervals)

plot(x, y) 
curve(coef1[pos] - best_k * coef2[pos] + coef2[pos] * x,
      from=min(x), to=best_k, add = TRUE)
segments(best_k, coef1[pos], max(x), coef1[pos])

- realcoolclub

---

你的数据长什么样？ - nrussell

@nrussell x = c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1)
y = c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516, 0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527, 0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.025414608) 我知道这些数据并不适合我所描述的拟合，但这也是我的重点... - realcoolclub

Answer 1

在Stack Overflow上有一个非常相似的帖子：Piecewise regression with a quadratic polynomial and a straight line joining smoothly at a break point。唯一的区别是现在我们考虑：

事实证明，在我的答案中定义的函数est、choose.c和pred根本不需要改变；我们只需要修改getX来返回您分段回归的设计矩阵：

getX <- function (x, c) cbind("beta0" = 1, "beta1" = pmin(x - c, 0))

现在，我们按照玩具示例中的代码来拟合您的数据模型：

x <- c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1)
y <- c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516, 
0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527, 
0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.25414608)

x 的取值范围为1到6，因此我们考虑

c.grid <- seq(1.1, 5.9, 0.05)
fit <- choose.c(x, y, c.grid)
fit$c
# 4.5

最后，我们制作预测图：

x.new <- seq(1, 6, by = 0.1)
p <- pred(fit, x.new)
plot(x, y, ylim = c(0, 0.4))
matlines(x.new, p[,-2], col = c(1,2,2), lty = c(1,2,2), lwd = 2)

我们在拟合模型中获得了丰富的信息:

str(fit)
#List of 12
# $ coefficients : num [1:2] 0.304 0.055
# $ residuals    : num [1:15] -0.06981 -0.08307 -0.02844 -0.00731 0.00624 ...
# $ fitted.values: num [1:15] 0.112 0.167 0.222 0.112 0.167 ...
# $ R            : num [1:2, 1:2] -3.873 0.258 9.295 -4.37
# $ sig2         : num 0.00401
# $ coef.table   : num [1:2, 1:4] 0.3041 0.055 0.0384 0.0145 7.917 ...
#  ..- attr(*, "dimnames")=List of 2
#  .. ..$ : chr [1:2] "beta0" "beta1"
#  .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
# $ aic          : num -34.2
# $ bic          : num -39.5
# $ c            : num 4.5
# $ RSS          : num 0.0521
# $ r.squared    : num 0.526
# $ adj.r.squared: num 0.49

例如，我们可以检查系数摘要表：

fit$coef.table
#        Estimate Std. Error  t value     Pr(>|t|)
#beta0 0.30406634 0.03840657 7.917039 2.506043e-06
#beta1 0.05500095 0.01448188 3.797915 2.216095e-03

- Zheyuan Li

Answer 2

今日免费次数已满, 请开通会员/明日再来

x2 = (x-k)*(x>k)
lm( y ~ x2)

或者，您可以使用I()

。

lm(y~ I((x-k)*(x>k)))

I()会逐字接受其中的内容，并忽略其他可能存在的（误）解释和其所在函数的功能。

如果您没有明确定义的k，那么您将需要在不同的k值上优化类似离差的东西。

- Max Candocia

---

我对如何将半条直线的斜率设置为0感到困惑？ - realcoolclub

一半的线总是为0，因为当(x<=k) 时，(x>k) 为FALSE，且 FALSE * anynumber 是0。当 x==k 时，它仍然为零，因为函数应该是连续的。之后，(x-k) 将以恒定速率增加。 - Max Candocia

啊，原来是这样！谢谢！ - realcoolclub