分段包断点是可变的，并在断点上查找标准错误答案

【问题标题】：Segmented package breakpoints are variable and finding standard errors on breakpoints分段包断点是可变的，并在断点上查找标准错误
【发布时间】：2020-06-10 15:04:36
【问题描述】：

今天有两个问题：一个是关于使用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]]

【问题讨论】：

一些可能的澄清：（1）gbay$pop 没有被使用，所以你可以删除它，（2）记住library(ggplot2) 和library(segmented)，（3）不需要做glm(..., family = gaussian())。这与lm(...). 相同，(4) 也许清楚你对“y 中的错误”的意思是什么？变化点是 x 轴上的一个位置。是预测 y 的不确定性吗？

标签： r linear-regression piecewise

【解决方案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的作者。

【讨论】：