分段程序包断点是可变的，并且在断点上发现标准错误

Question

[今天有两个问题：一个问题是关于使用segmented包创建分段回归并在我多次运行模型时获得不同的断点，第二个问题是关于断点的标准误差。

加载并查看数据：

gbay<-data.frame(matrix(,nrow=46,ncol=3))
colnames(gbay)<-c("pop","cal.length","temp")

gbay$cal.length<-c(0.597, 0.834, 1.120, 1.353, 0.119, 1.301, 0.944, 3.127, 3.375, 3.171, 3.400, 3.376, 3.322, 3.785, 3.304, 3.197, 3.216,
 4.183, 2.171, 3.989, 3.187 ,4.153, 3.252, 4.960, 4.268, 4.827, 4.869, 3.932, 4.573, 4.645, 4.634, 4.713, 4.879, 4.724,
5.031, 4.746, 5.047, 5.714, 5.227, 4.701,4.280, 5.296, 4.977, 4.932, 4.374, 4.758)

gbay$temp<-c(16, 16, 16, 16, 16, 16, 16, 20, 20, 20, 20, 20, 20, 20, 20, 24, 24, 24, 24, 24, 24, 24, 24, 26, 26, 26, 26, 26, 26, 26, 28, 28, 28, 28,
28, 28, 28, 28, 28, 30, 30, 30, 30, 30, 30, 30)
gbay$pop<-gb

ggplot(gbay,aes(x=temp,y=cal.length))+geom_point()

1）以上是我在实验室温度范围内附加的蜗牛生长速率数据。我想创建一个分段回归，该回归将确定具有断点的最佳生长温度和最大生长速率。我一直使用segmented包取得了一些成功，但是，我的数据遇到了一些困难，因为该包在两个断点之间来回移动，一个在27.56（根据原始数据和我的问题是我想要的断点x组件），另一个是在18.59，发生在我不希望模型计算的“假”最优条件下。我曾尝试在分段模型中操纵psi，但每次运行模型时都会得到50-50的断点。有没有办法告诉程序包只关注某些x边界来搜索断点？

m.gbay<-glm(cal.length~temp,gbay,family=gaussian)
seg.gbay<-segmented(m.gbay,seg.Z = ~temp, psi=28)
xmin<-min(gbay$temp,na.rm=T)
xmax<-max(gbay$temp,na.rm=T)
predicted.gbay<-data.frame(temp=seq(xmin,xmax,length.out=100))
predicted.gbay$cal.length<-predict(seg.gbay,predicted.gbay)
predicted.gbay$pop<-"gb"

ggplot(predicted.gbay,aes(x=temp,y=cal.length))+geom_line(aes(x=temp,y=cal.length))+
  ylab("Shell Length (mm)")+xlab("Common Garden Temperature (°C)")

summary(seg.gbay)

2）我正在尝试从此数据中提取断点（psi）的x和y分量。我已经成功完成了。但是，我也希望能够为断点提取x和y分量的误差。我认为模型会吐出x分量的标准错误，但我想知道segmented或其他程序包中是否有办法在断点的y分量中找到错误？

breakpts<-data.frame(matrix(,nrow=1,ncol=4))
colnames(breakpts)<-c("brkptx","brkpty","x_err","y_err")

breakpts[1,1]<-seg.gbay$psi[[2]]
breakpts[1,2]<-(seg.gbay$psi[[2]]*coef(seg.gbay)[[2]])+(coef(seg.gbay)[[1]])
breakpts[1,3]<-seg.gbay$psi[[3]]

Answer 1

此数据集非常小，因此，如果您没有先验知识来指导它，那么更改点将位于不明确的位置。大多数现有软件包仅标识一个更改点，而没有量化它们通常很奇怪的分布。

我认为mcp软件包可以满足您的需求。简要地说，您先对一条直线进行建模，然后再对其进行连接：

model = list(
  cal.length ~ 1 + temp,  # Line with intercept
  ~ 0 + temp  # Joined slope
)

现在让我们适应它。 family = gaussian()是隐式的。请注意，我设置了很多迭代和并行处理，因为在这种情况下更改点的位置非常难以识别，因此需要大量工作来探索后验：

library(mcp)
fit = mcp(model, data = gbay, iter = 100000, cores = 3)

默认图显示了估计变化点的后验（蓝线）。我们可以添加拟合间隔（红色）和预测间隔（绿色）。灰线是来自后验的25个随机样本：

plot(fit, q_fit = T, q_predict = T)

如您所见，变更点的后验是双峰的。您也可以调用plot_pars(fit)和summary(fit)查看有关各个参数的更多详细信息。如果要测试最早模式与第二种模式的证据，可以使用hypothesis(fit, "cp_1 < 25")。如果您具有有关可靠坡度，更改点位置等的先验知识，则可以轻松地使用例如mcp(model, gbay, prior = list(cp_1 = "dunif(0, 25)"))进行添加。

在mcp website和this preprint中了解有关mcp的更多信息，包括安装说明。免责声明：我是mcp的作者。