我有 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“测试”。 ...但也许还有其他想法?
谢谢!非常感激!