如何拟合分段常数(或阶跃函数)模型并与 R 中的逻辑模型进行比较

机器算法验证 回归 物流 模型 aic nls
2022-03-06 16:28:56

我有 x,y 数据,其中 x 是位置(沿样带),y 是连续变量(例如温度)。

x<-c(115,116,117,118,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151)
y<-c(68.3864344929207,69.3317985238957,70.5644430507252,68.7245413735345,68.2444929220624,69.2335442981487,68.9683874048675,69.7104166614154,70.6047681836632,71.1807115091657,70.3333333333333,70,69.735105818315,69.6487585375593,70,69.4943374527374,69,70.3590104484011,70.283981899468,68.9087146344906,68.6666666666667,68.5232663969658,68.088410889283,67.567883666713,66.9865494087022,66,66.3333333333333,66.0165138638852,65.6666666666667,65.7465975732667,66.3333333333333,66.6579079386334)

我想拟合一个线性模型、一个 4 参数逻辑模型和一个分段常数函数(又名阶跃函数),这些数据有一个断点,并比较模型以确定 y 是否逐渐变化,比线性变化更突然,还是一步突然变化-喜欢时尚。

(注意:绘制数据表明逻辑模型在这种情况下是最合适的,但我最终将不得不在没有这些信息的情况下提前遍历许多不同的因变量)。

我将 nls() 与手动公式 (y~int+slope) 用于线性模型,将 SSfpl 用于 4 参数逻辑。

l<-nls(y~int+x*slope,data=data.frame(cbind(x,y)),start=list(int=0,slope=1)) # linear model 
fl<-nls(y~SSfpl(x,A,B,xmid,scal),data=data.frame(cbind(x,y))) # four parameter logistic

我的问题是拟合分段常数模型。似乎存在用于分段线性模型的函​​数,但用于阶跃函数的函数并不多。我发现:

library(cumSeg)
pc<-jumpoints(y,k=1,output="1")

这给了我 x 中最佳断点的索引:

breakpt<-pc$psi

但我不确定如何将输出与线性或逻辑模型进行比较。例如,我认为我必须使用 AIC 进行模型比较,但不确定如何根据 jumppoints() 返回的对象的类来计算它。所以我手动尝试使用最佳断点生成模型:

pcm<-lm(y~I(x<=x[breakpt])+I(x>x[breakpt]))

现在我可以获得所有模型的 AIC 分数。但我注意到,当我绘制手动步进函数 (pcm) 与 cumSeg (pc) 中生成的拟合时,看起来我的手动模型的拟合效果不太好......

par(mfrow=c(1,2))
plot(x,y)
lines(x,predict(cpm) # versus
plot(pc)

想知道为什么会这样?

并寻找解决整体问题的替代建议。在模型拟合方面:主要问题是我如何判断一个变量是否是阶梯式的,而不是其他或多或少线性变化?(我可以有一个比线性模型更适合的逻辑模型,但在 xmid 有类似的斜率......所以这不是一个真正的答案)。

在模型比较结束时:我使用 MuMIn 包中的 AICc() 函数来获取不同模型的 AIC 值,并使用 qpcR 包中的 akaike.weights() 函数来执行增量 AIC“测试”。 ...但也许还有其他想法?

谢谢!非常感激!

1个回答

我认为主要困难之一是分段常数回归通常不称为“分段常数回归”。它们通常被称为回归树,这是一个很好的视觉名称,但如果您还不知道人们如何称呼它们,则不是特别适合用 google 搜索!它们可以R与内置的rpart包配合使用(我相信rpart代表“递归分区”,这是我们对同一概念的第三个名称)。

以下是rpart对您的数据的实际操作:

df <- data.frame(x=x, y=y)
tree <- rpart(y ~ x, data=df)

我写了一个plot_tree显示预测的小函数

plot_tree <- function(tree, x, y) {
  s <- seq(110, 155, by=.5)
  plot(x, y)
  lines(s, predict(tree, data.frame(x=s)))
}

当应用于数据的默认树时,它看起来像这样

plot_tree(tree, x, y)

在此处输入图像描述

您可以通过使用来控制适合数据的粒度rpart.controll

tree <- rpart(y ~ x, data=df, control=rpart.control(minsplit=5, cp=.0001))
plot_tree(tree, x, y)

在此处输入图像描述