ggplot 显示自举曲线拟合的置信区间答案

【问题标题】：ggplot to show confidence intervals from bootstrapping curve fittingggplot 显示自举曲线拟合的置信区间
【发布时间】：2017-03-20 00:01:05
【问题描述】：

感谢here 的帮助，我能够使用自举法获得曲线拟合的意大利面条图。我试图从这些拟合模型中得出置信带。我没有运气得到类似的东西

quants <- apply(fitted_boot, 1, quantile, c(0.025, 0.5, 0.975))

使用以下内容：

library(dplyr)
library(broom)
library(ggplot2)

xdata <- c(-35.98, -34.74, -33.46, -32.04, -30.86, -29.64, -28.50, -27.29, -26.00, 
           -24.77, -23.57, -22.21, -21.19, -20.16, -18.77, -17.57, -16.47, -15.35,
           -14.40, -13.09, -11.90, -10.47, -9.95,-8.90,-7.77,-6.80, -5.99,
           -5.17, -4.21, -3.06, -2.29, -1.04)
ydata <- c(-4.425, -4.134, -5.145, -5.411, -6.711, -7.725, -8.087, -9.059, -10.657,
           -11.734, NA, -12.803, -12.906, -12.460, -12.128, -11.667, -10.947, -10.294,
           -9.185, -8.620, -8.025, -7.493, -6.713, -6.503, -6.316, -5.662, -5.734, -4.984,
           -4.723, -4.753, -4.503, -4.200)

data <- data.frame(xdata,ydata)
x_range <- seq(min(xdata), max(xdata), length.out = 1000)

fitted_boot <- data %>% 
  bootstrap(100) %>%
  do({
    m <- nls(ydata ~ A*cos(2*pi*((xdata-x_0)/z))+M, ., start=list(A=4,M=-7,x_0=-10,z=30))
    f <- predict(m, newdata = list(xdata = x_range))
    data.frame(xdata = x_range, .fitted = f)
    } )

ggplot(data, aes(xdata, ydata)) +
  geom_line(aes(y=.fitted, group=replicate), fitted_boot, alpha=.1, color="blue") +
  geom_point(size=3) +
  theme_bw()

我认为 geom_ribbon() 可能是一个不错的方法，但我只是不知道从这里去哪里。

感谢 Axeman 对其他帖子的帮助！

【问题讨论】：

标签： r ggplot2 confidence-interval

【解决方案1】：

一种方法是计算每个 x 值的置信区间，然后将其绘制出来。在这里，我使用了第 2.5 个百分位和第 97.5 个百分位之外的第一个值，尽管您可以根据需要调整代码。

首先，我将group_by 更改为xdata 位置（而不是复制）。然后，我通过.fitted 值arrange，以便我可以slice 取出我想要的值（百分位截止值之外的第一个值）。最后，我用我得到的界限标记它们（它们总是先低于后高于，因为我们排序了）。

forConfInt <-
  fitted_boot %>%
  ungroup() %>%
  group_by(xdata) %>%
  arrange(.fitted) %>%
  slice(c(floor(0.025 * n() )
          , ceiling(0.975 * n() ) ) ) %>%
  mutate(range = c("lower", "upper"))

这给出了：

   replicate     xdata   .fitted range
       <int>     <dbl>     <dbl> <chr>
1          9 -35.98000 -4.927462 lower
2         94 -35.98000 -4.249348 upper
3          9 -35.94503 -4.927248 lower
4         94 -35.94503 -4.257776 upper
5          9 -35.91005 -4.927228 lower
6         94 -35.91005 -4.266334 upper
7          9 -35.87508 -4.927401 lower
8         94 -35.87508 -4.275020 upper
9          9 -35.84010 -4.927766 lower
10        94 -35.84010 -4.283836 upper
# ... with 1,990 more rows

然后我们可以在ggplot 调用中添加一行：

ggplot(data, aes(xdata, ydata)) +
  geom_line(aes(y=.fitted, group=replicate), fitted_boot, alpha=.1, color="blue") +
  # Added confidence interval:
  geom_line(aes(y=.fitted, group=range), forConfInt, color="red") +
  geom_point(size=3) +
  theme_bw()

给出这个情节：

【讨论】：