在 R 中绘制具有 95% 置信区间的阈值/分段/变化点模型答案

【问题标题】：plotting threshold/piecewise/change point models with 95% confidence intervals in R在 R 中绘制具有 95% 置信区间的阈值/分段/变化点模型
【发布时间】：2020-05-29 21:54:52
【问题描述】：

我想绘制一个阈值模型，线段之间具有平滑的 95% 置信区间线。您会认为这很简单，但我一直无法找到答案！

我的阈值/断点是已知的，如果有办法可视化这些数据，那就太好了。我已经尝试了产生以下情节的分段包：

该图显示了一个阈值模型，其断点位于 5.4。但是，回归线之间的置信区间并不平滑。

如果有人知道在分段回归线（理想情况下在 ggplot 中）之间产生平滑（即没有线段之间的跳跃）CI 线的任何方法，那将是惊人的。非常感谢。

我已经包含了示例数据和我在下面尝试过的代码：

x <- c(2.26, 1.95, 1.59, 1.81, 2.01, 1.63, 1.62, 1.19, 1.41, 1.35, 1.32, 1.52, 1.10, 1.12, 1.11, 1.14, 1.23, 1.05, 0.95, 1.30, 0.79,
0.81, 1.15, 1.10, 1.29, 0.97, 1.05, 1.05, 0.84, 0.64, 0.80, 0.81, 0.61, 0.71, 0.75, 0.30, 0.30, 0.49, 1.13, 0.55, 0.77, 0.51,
0.67, 0.43, 1.11, 0.29, 0.36, 0.57, 0.02, 0.22, 3.18, 3.79, 2.49, 2.44, 2.12, 2.45, 3.22, 3.44, 3.86, 3.53, 3.13)

y <- c(22.37, 18.93, 16.99, 15.65, 14.62, 13.79, 13.09, 12.49, 11.95, 11.48, 11.05, 10.66, 10.30,  9.96,  9.65,  9.35,  9.07,  8.81,
       8.56,  8.32,  8.09,  7.87,  7.65,  7.45,  7.25,  7.05,  6.86,  6.68,  6.50,  6.32,  6.15,  5.97,  5.80,  5.63,  5.47,  5.30,
        5.13,  4.96,  4.80,  4.63,  4.45,  4.28,  4.09,  3.90,  3.71,  3.50,  3.27,  3.01,  2.70,  2.28, 22.37, 16.99, 11.05,  8.81,
       8.56,  8.32,  7.25,  7.05,  6.50,  6.15,  5.63)

lin.mod <- lm(y ~  x)
segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=2)
plot(x, y)
plot(segmented.mod, add=TRUE, conf.level = 0.95)

产生以下图（以及 95% 置信区间的相关跳跃）：

segmented plot

【问题讨论】：

问题在于，在分段回归中，置信区间不是平滑和连续的。如果你想显示 95% 的置信区间，这就是它们的样子。如果您想显示平滑的线条，那么绘图可能看起来更好，但您将显示的不是 95% 的置信区间。如果您需要帮助使线条平滑，那当然可以，但我们需要查看一些示例数据以及您迄今为止所做的尝试才能为您提供帮助。
@AllanCameron 非常感谢，这当然很有意义。我现在已经包含了一些示例数据以及生成的图。我很想听听关于制作平滑线的任何想法（当然要承认这些不再是 95% 的置信区间）。再次感谢您，非常感谢您的帮助！
@AllanCameron，为什么 CI 不应该在更改点周围平滑？看我的回答。
@JonasLindeløv 我真的很喜欢你的回答（赞成），我认为固定截止的第二部分最接近 OP 所寻找的。我认为您的第一种方法总体上可能是一种更好的统计方法，尽管如果没有对问题的清晰描述就很难知道。过去，当我想做分段回归时，我的方法是找到一个非线性模型。这有时会导致更好地了解手头的问题 - 例如参见dx.doi.org/10.1136/annrheumdis-2013-203293。但我认为你的回答应该在这里被接受。

标签： r ggplot2 threshold piecewise

【解决方案1】：

背景：现有变更点包中的不平滑是由于常客包以固定的变更点值运行。但与所有推断参数一样，这是错误的，因为变化的位置确实存在不确定性。

解决方案： AFAIK，只有贝叶斯方法可以量化这一点，mcp 包填补了这个空间。

library(mcp)
model = list(
  y ~ 1 + x,   # Segment 1: Intercept and slope
  ~ 0 + x  # Segment 2: Joined slope (no intercept change)
)
fit = mcp(model, data = data.frame(x, y))

默认绘图（plot.mcpfit() 返回 ggplot 对象）：

plot(fit) + ggtitle("Default plot")

每一行代表生成数据的可能模型。变化点的后验显示为蓝色密度。您可以使用plot(fit, q_fit = TRUE) 在顶部添加一个可信区间或单独绘制它：

plot(fit, lines = 0, q_fit = c(0.025, 0.975), cp_dens = FALSE) + ggtitle("Credible interval only")

如果您的更改点确实已知，并且如果您想为每个段建模不同的残差比例（即，准模拟segmented），您可以这样做：

model2 = list(
  y ~ 1 + x,
  ~ 0 + x + sigma(1)  # Add intercept change in residual scale
)
fit = mcp(model2, df, prior = list(cp_1 = 1.9))  # Note: prior is a fixed value - not a distribution.
plot(fit, q_fit = TRUE, cp_dens = FALSE)

请注意，CI 不会像 segmented 那样在更改点周围“跳跃”。我相信这是正确的行为。披露：我是mcp的作者。

【讨论】：

非常感谢！这非常有帮助，这正是我所追求的。在我的情况下，通过比较不同的阈值模型并将顶级模型（AIC）作为变化点，变化点是“已知的”。但是现在我看到 mcp 包可以做什么，我更倾向于使用您在此处提出的第一种方法。再次感谢！
再问几个关于 mcp 包的问题，以防你碰巧再次访问这篇文章。有没有办法在代表每个可能模型的线中取一条平均线？这让我想到了我的第二个问题:) 是否可以访问原始拟合模型数据来自定义绘图？ plot.mcpfit 和 plot_pars 很棒，只是想知道是否可以访问原始数据（如果我错过了这一点，很抱歉，这很明显）。谢谢！
@Emmax，您可以使用plot(fit, q_fit = c(0.5)) 绘制中位数（50% 分位数）。原始样本以mcmc.list 的形式存储在fit$samples 中，这是一种标准格式，可以由许多包处理。 mcp 经常使用 tidybayes 和 bayesplot。 See more here.