R中的分段回归：绘制线段

Question

我有54个点。它们代表产品的供需情况。我想要展示在供应方面存在一个突破点。

首先，我对x轴（供应）进行排序并删除出现两次的值。我有47个值，但我删除了第一个和最后一个（考虑它们作为突破点没有意义）。突破长度为45：

Break<-(sort(unique(offer))[2:46])

然后，对于这些可能的断点，我估计一个模型，并将“d”中的残差标准误差（模型摘要对象中的第六个元素）保留下来。

d<-numeric(45)
for (i in 1:45) {
model<-lm(demand~(offer<Break[i])*offer + (offer>=Break[i])*offer)
d[i]<-summary(model)[[6]] }

绘制d后，我注意到我的较小残差标准误差是34，对应于"Break[34]": 22.4。因此我用我的最终断点编写了模型：

model<-lm(demand~(offer<22.4)*offer + (offer>=22.4)*offer)

最后，我对我的新模型感到满意。它比简单的线性模型显着优秀。现在我想画出来：

plot(demand~offer)
i <- order(offer)
lines(offer[i], predict(model,list(offer))[i])

但是我收到一个警告信息：

Warning message:
In predict.lm(model, list(offer)) :
  prediction from a rank-deficient fit may be misleading

而且更重要的是，我的图中的线条非常奇怪。

这是我的数据：

demand <- c(1155, 362, 357, 111, 703, 494, 410, 63, 616, 468, 973, 235,
            180, 69, 305, 106, 155, 422, 44, 1008, 225, 321, 1001, 531, 143,
            251, 216, 57, 146, 226, 169, 32, 75, 102, 4, 68, 102, 462, 295,
            196, 50, 739, 287, 226, 706, 127, 85, 234, 153, 4, 373, 54, 81,
            18)
offer <- c(39.3, 23.5, 22.4, 6.1, 35.9, 35.5, 23.2, 9.1, 27.5, 28.6, 41.3,
           16.9, 18.2, 9, 28.6, 12.7, 11.8, 27.9, 21.6, 45.9, 11.4, 16.6,
           40.7, 22.4, 17.4, 14.3, 14.6, 6.6, 10.6, 14.3, 3.4, 5.1, 4.1,
           4.1, 1.7, 7.5, 7.8, 22.6, 8.6, 7.7, 7.8, 34.7, 15.6, 18.5, 35,
           16.5, 11.3, 7.7, 14.8, 2, 12.4, 9.2, 11.8, 3.9)

- Antonin

---

54个点似乎不是一个很大的点数来检测这样的转变。您可以选择在stats.stackexchange.com上发布此问题，并具体说明是否有足够的点来检测数据中的突变。这只是我的两分钱。 - Paul Hiemstra

我认为这在统计上相当可疑。最好在模型本身中估计断点（尽管这使其非线性）。您不能信任当前非正式估计过程中的p-值或标准误差。 - hadley

54分并不算很多，我同意，但是我的线性回归和分段线性回归都是显著的。而且，与自由度少两个的线性模型相比，分段线性模型的残差标准误差为103.9，而线性模型为121.3。分段模型显著更好。 - Antonin

Answer 1

这里有一种更简单的方法，使用ggplot2。

require(ggplot2)
qplot(offer, demand, group = offer > 22.4, geom = c('point', 'smooth'), 
   method = 'lm', se = F, data = dat)

编辑：我还建议查看这个包segmented，它支持自动检测和估计分段回归模型。

更新：

以下是使用R包segmented自动检测断点的示例。

library(segmented)
set.seed(12)
xx <- 1:100
zz <- runif(100)
yy <- 2 + 1.5*pmax(xx - 35, 0) - 1.5*pmax(xx - 70, 0) + 15*pmax(zz - .5, 0) + 
  rnorm(100,0,2)
dati <- data.frame(x = xx, y = yy, z = zz)
out.lm <- lm(y ~ x, data = dati)
o <- segmented(out.lm, seg.Z = ~x, psi = list(x = c(30,60)),
  control = seg.control(display = FALSE)
)
dat2 = data.frame(x = xx, y = broken.line(o)$fit)

library(ggplot2)
ggplot(dati, aes(x = x, y = y)) +
  geom_point() +
  geom_line(data = dat2, color = 'blue')

- Ramnath

---

感谢您提出使用“segmented”软件包的想法。“Muggeo，V.M.R.（2003）Estimating regression models with unknown break-points. Statistics in Medicine 22, 3055–3071”是一篇有趣的论文，可以帮助理解该软件包中发生的情况。 - Antonin

特别是，它对我使用的代码有一个优势：这两个线段是连接的！在《R语言实战》一书中，作者没有提到他的线段不是连接的，甚至展示了一个带有连接线段的图表... - Antonin

使用ggplot()和segmented()的示例在哪里？我似乎无论在哪里都找不到。 - Adam Erickson

我已添加了一个使用“ggplot”和“segmented”的示例。 - Ramnath

在使用 qplot(…) 后，我得到了一个错误：Error: Unknown parameters: method, se。 - sam

你的解决方案也适用于这里吗？ - rnorouzian

Answer 2

已经为你指明了正确的方向。在你得到的图形中，唯一“奇怪”的是lines在每个连续点之间画出一条直线，这意味着你看到的“跳跃”只是连接每条线段两端的结果。

如果你不想要这种连接符号，你需要将lines拆分成两个独立的部分。

此外，我认为你可以简化你的回归曲线。以下是我的方法:

#After reading your data into dat
Break <- 22.4
dat$grp <- dat$offer < Break

#Note the addition of the grp variable makes this a bit easier to read
m <- lm(demand~offer*grp,data = dat)
dat$pred <- predict(m)

plot(dat$offer,dat$demand)
dat <- dat[order(dat$offer),]
with(subset(dat,offer < Break),lines(offer,pred))
with(subset(dat,offer >= Break),lines(offer,pred))

这将产生以下绘图：

- joran

Answer 3

奇怪的线条只是由于点绘制的顺序不同造成的。 下面的绘图应该会更好看：

i <- order(offer)
lines(offer[i], predict(model,list(offer))[i])

警告是因为*字符被lm解释了。

> lm(demand~(offer<22.4)*offer + (offer>=22.4)*offer)
Call:
lm(formula = demand ~ (offer < 22.4) * offer + (offer >= 22.4) * offer)
Coefficients:
            (Intercept)         offer < 22.4TRUE                    offer  
                -309.46                   356.08                    29.86  
      offer >= 22.4TRUE   offer < 22.4TRUE:offer  offer:offer >= 22.4TRUE  
                     NA                   -20.79                       NA  

此外，(offer<22.4)*offer是一个不连续的函数：这就是不连续性的来源。

以下内容应更符合您的要求。

model <- lm(
  demand ~ ifelse(offer<22.4,offer-22.4,0) 
           + ifelse(offer>=22.4,offer-22.4,0) )

- Vincent Zoonekynd

---

感谢回答！我仍然收到警告信息：“警告信息：在predict.lm(model2, list(offer))中： 预测来自秩缺陷拟合可能是误导性的。” 结果稍微好一些了：我只有一个分段，但我不明白为什么会有3个分段（即，我的两个分段不是自然连接的...） - Antonin

请问您能否解释一下为什么这会使它更平滑？直觉上认为应该是函数的第一部分是22.4-offer，第二部分是offer-22.4。但是第一部分不是意味着“如果报价低于22.4，则计算出的报价减去22.4”，这意味着它将给出负的协变量吗？因此，当我们估计输入时，我们的系数将乘以这些负值吗？ - chris