分段包断点是可变的，在断点上查找标准错误
Segmented package breakpoints are variable and finding standard errors on breakpoints
今天有两个问题：一个是关于使用 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，它发生在我不希望模型计算的 "false" 最优值。我试过在分段模型中操纵 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]]
这个数据集非常小，因此如果您没有先验知识来指导它，那么变化点所在的位置是不明确的。大多数现有的包只是识别一个单一的变化点，而没有量化它们通常很奇怪的分布。
我认为 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 and in this preprint 上阅读有关 mcp 的更多信息，包括安装说明。披露：我是 mcp.
的作者
今天有两个问题：一个是关于使用 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，它发生在我不希望模型计算的 "false" 最优值。我试过在分段模型中操纵 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]]
这个数据集非常小，因此如果您没有先验知识来指导它，那么变化点所在的位置是不明确的。大多数现有的包只是识别一个单一的变化点，而没有量化它们通常很奇怪的分布。
我认为 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 and in this preprint 上阅读有关 mcp 的更多信息，包括安装说明。披露：我是 mcp.