难以在 R 中拟合分段线性数据

Posted 2023-02-14

【中文标题】难以在 R 中拟合分段线性数据

【英文标题】：Difficulty fitting piecewise linear data in R

【发布时间】：2022-01-05 13:06:51

【问题描述】：

我有以下数据（产品成本与时间），如下所示：

annum <- c(1903, 1904, 1905, 1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 
    1914, 1915, 1916, 1917, 1918, 1919)
cost <- c(0.0000,  18.6140,  92.1278, 101.9393, 112.0808, 122.5521, 
    133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
    339.6527, 362.3537, 377.7775, 402.8443, 437.5539)

mydata <- as.data.frame(cbind(annum, cost))

g <- ggplot(mydata, aes(x = annum, y = cost))
g <- g + geom_point()
g <- g + scale_y_continuous(labels=scales::dollar_format())
g

This is the resulting plot of this data using this code 该图显示了对我来说看起来是分段线性的东西；从 1904 年到 1905 年有一个台阶；然后是从 1905 年到 1910 年的清晰界限；然后一步；然后是从 1911 到结尾的另一行。 （第一点 (1903, 0) 是虚构的。）

我尝试使用分段包对此进行建模，但它没有选择像 1904.5 和 1910.5 这样的断点，而是在 1911 和 1912 之间找到两个点。

我尝试了一些其他技术（例如，“The R Book”中的“蛮力”和直接拟合），但我显然没有达到我需要的程度。任何帮助将不胜感激。

理想情况下，我最终会得到每个段的方程和显示分段拟合和拟合置信区间的单个图。

【参考方案1】：

可以为此使用包 struccchange。这里有一个简化的代码版本：

library("strucchange")

startyear <- startyear
cost <- c(0.0000,  18.6140,  92.1278, 101.9393, 112.0808, 122.5521, 
          133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
          339.6527, 362.3537, 377.7775, 402.8443, 437.5539)

ts <- ts(cost, start=1903)
plot(ts)

## for small data sets you might consider to reduce segment length
bp <- breakpoints(ts ~ time(ts), data=ts, h = 5)

## BIC selection of breakpoints
plot(bp)
breakdates(bp)
fm1 <- lm(ts ~ time(ts) * breakfactor(bp), data=ts)
coef(fm1)

plot(ts, type="p")
lines(ts(fitted(fm1),  start = startyear),  col = 4)
lines(bp)
confint(bp)

lines(confint(bp))

更多信息可以在包装小插图或相关出版物之一中找到，例如https://doi.org/10.18637/jss.v007.i02 因此，例如可以进行显着性检验、估计置信区间或包含协变量。

段长度为 2 是不可能的，因为无法估计残差。同样，只有当段足够长时，才能估计置信区间。因此，下面只显示了一个断点，而@Rui Barradas 的优秀答案省略了置信区间但显示了两个断点。

她是一个没有前两点的例子和一个额外的假设来估计小段的置信区间：

library("strucchange")

startyear <- 1905
cost <- c(92.1278, 101.9393, 112.0808, 122.5521, 
          133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
          339.6527, 362.3537, 377.7775, 402.8443, 437.5539)

ts <- ts(cost, start=startyear)
bp <- breakpoints(ts ~ time(ts), data=ts, h = 5)
fm1 <- lm(ts ~ time(ts) * breakfactor(bp), data=ts)
plot(ts, type="p")
lines(ts(fitted(fm1),  start = startyear),  col = 4)
lines(confint(bp, het.err=FALSE))

编辑：

修正了原始版本的错误 添加了系数和置信区间 已添加图片 添加了省略的前 2 个值的示例

【讨论】：

当我删除前两个点（第一个是虚构的）时，早期段的拟合失败（斜率不正确）。对此有什么想法吗？

感谢您的评论。原始版本包含 2 个错误。最重要的是在lm模型公式中使用*，另一个是时间序列对象和原始向量的混合。

【参考方案2】：

这是另一个使用包strucchange 的解决方案，但没有先创建时间序列。

library(strucchange)

# first get a segment size as a fraction 
# of the number of observations
n <- nrow(mydata)
segmts <- 3
h <- (segmts + 1)/n

# now estimate the breakpoints
b <- breakpoints(cost ~ annum, h = h, breaks = (segmts - 1L), data = mydata)
bp <- mydata[b$breakpoints, "annum"]

# create a grouping variable for `ggplot`
# each group is a segment
bp <- c(bp, Inf)
mydata$grp <- findInterval(mydata$annum, bp, left.open = TRUE)

# plot the linear regressions
g + geom_smooth(
  mapping = aes(group = grp),
  method = "lm",
  formula = y ~ x,
  se = FALSE
)

---

如果第一个数据点被删除，将只有两个段，但上面的代码仍然可以工作。

mydata <- mydata[-(1:2), ]
n <- nrow(mydata)
segmts <- 2
h <- (segmts + 1)/n
b <- breakpoints(cost ~ annum, h = h, breaks = segmts - 1L, data = mydata)
bp <- mydata[b$breakpoints, "annum"]
bp <- c(bp, Inf)
mydata$grp <- findInterval(mydata$annum, bp, left.open = TRUE)
mydata$grp <- factor(mydata$grp)

g + geom_smooth(
  mapping = aes(group = grp),
  method = "lm",
  formula = y ~ x,
  se = FALSE
)

【讨论】：

如果我切断前两点（无论如何，第一点都是虚构的），这只会留下两段。当我运行这个方法时，设置segmnts <- 2，我得到错误：Error in breakpoints.formula(Total ~ Year, h = h, breaks = (segmts - 1L), : minimum segment size must be greater than the number of regressors.有什么想法吗？

试试h <- (segmts + 1)/n。它在没有前 2 点的情况下有效。我将编辑我的答案。

【参考方案3】：

变化点问题的置信区间对于频率论方法来说是一个难题，例如strucchange。通常，您只需获得每个段的置信区间，即段之间的硬中断而不是平滑过渡。

使用贝叶斯方法更简单。这是使用mcp 包的解决方案。为了炫耀，我们绘制了拟合区间和（红色虚线）和预测区间（绿色虚线）。灰线是从后验分布中随机抽取的，x 轴上的密度是变化点位置的后验。

data = data.frame(
  annum = 1903:1919,
  cost = c(0.0000,  18.6140,  92.1278, 101.9393, 112.0808, 122.5521, 
          133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
          339.6527, 362.3537, 377.7775, 402.8443, 437.5539)
)

# Model as three disjoined slopes
model = list(
  cost ~ 1 + annum,
  ~ 1 + annum,
  ~ 1 + annum
)

library(mcp)
fit = mcp(model, data)
plot(fit, q_fit = TRUE, q_predict = TRUE)

如果您对变化点和段的参数估计感兴趣，请致电summary(fit)：

        name    mean  lower    upper Rhat n.eff
     annum_1   -0.11   -0.2 -6.6e-04  2.5    25
     annum_2   10.36    7.4  1.3e+01  1.0   609
     annum_3   22.74   21.2  2.4e+01  1.0   264
        cp_1 1904.50 1904.0  1.9e+03  2.5    24
        cp_2 1910.46 1910.0  1.9e+03  1.0   778
 Intercept_1  221.39   10.8  3.9e+02  1.0   948
 Intercept_2   86.77   75.0  9.8e+01  1.0  1297
 Intercept_3  236.03  221.7  2.5e+02  1.0   237
     sigma_1    5.97    3.6  8.9e+00  1.0  1709

【讨论】：

这看起来是一种有趣的方法，但无论出于何种原因，我都无法重现您所做的事情。我收到一个巨大的 JAGS 错误，开头为：`Error: .onLoad failed in loadNamespace() for 'rjags', details: call: dyn.load(file, DLLpath = DLLpath, ...) 错误：无法加载共享对象'/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rjags/libs/rjags.so'：我去了lindeloev网站尝试修复它，但没有任何乐趣。有什么想法吗？

您安装了 JAGS 吗？ sourceforge.net/projects/mcmc-jags/files/JAGS/4.x 谷歌搜索你的错误信息，这是其他人的问题：gist.github.com/casallas/8411082。 JAGS 是用于 MCMC 采样的软件。

太棒了！我知道这是题外话，但是您将如何使用 ggplot 进行绘制？

情节是ggplot，所以可以plot(fit) + labs(title = "This is the title")。要从头开始，这里有一个例子：lindeloev.github.io/mcp/articles/…。基本上，只需使用fitted(fit) 或fitted(fit, summary = FALSE) 作为数据，我相信你会弄明白的:-)

【参考方案4】：

这有帮助吗？使用loess方法？

library(tidyverse)
ggplot(mydata, aes(x = annum, y = cost))+
  geom_point()+
  geom_smooth(method = "loess", formula = "y~x")