如何将stat_ecdf与geom_ribbon结合起来？

Question

我正在尝试使用 ggplot2 绘制一些数据的 ECDF，其中“置信区间”通过阴影区域表示。我无法将

geom_ribbon()

与

stat_ecdf()

组合起来以达到我想要的效果。

考虑以下示例数据：

set.seed(1)
dat <- data.frame(variable = rlnorm(100) + 2)
dat <- transform(dat, lower = variable - 2, upper = variable + 2)

> head(dat)
  variable     lower    upper
1 2.534484 0.5344838 4.534484
2 3.201587 1.2015872 5.201587
3 2.433602 0.4336018 4.433602
4 6.929713 4.9297132 8.929713
5 3.390284 1.3902836 5.390284
6 2.440225 0.4402254 4.440225

我能够使用

生成

variable

的ECDF

library("ggplot2")
ggplot(dat, aes(x = variable)) +
    geom_step(stat = "ecdf")

但是，我无法使用

lower

和

upper

作为

ymin

和

ymax

的美学，将置信区间叠加在绘图上作为另一层。我试过了：geom_ribbon()

但这会引发以下错误

ggplot(dat, aes(x = variable)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), stat = "ecdf") +
    geom_step(stat = "ecdf")

有没有办法诱导 

Error: geom_ribbon requires the following missing aesthetics: ymin, ymax

与

geom_ribbon()

一起工作以产生阴影置信区间？或者，任何人都可以建议一种替代方法，将由

stat_ecdf()

和

lower

定义的阴影多边形作为 ECDF 图的图层添加到 ECDF 图中？

Answer 1

upper

好吧，这与您尝试做的事情不同，但它应该解释发生了什么。 

ggplot(dat, aes(x = variable)) +
  geom_ribbon(aes(x = variable,ymin = ..y..-2,ymax = ..y..+2), stat = "ecdf",alpha=0.2) +
  geom_step(stat = "ecdf")

返回一个数据框，其中仅包含原始 x 和计算出的 y，所以我认为这就是您需要处理的全部内容。即

stat

一次仅计算单个 x 的累积分布函数。我能想到的唯一的其他事情是显而易见的，分别计算下限和上限，如下所示：

stat_ecdf

Answer 2

l <- ecdf(dat$lower)
u <- ecdf(dat$upper)
v <- ecdf(dat$variable)
dat$lower1 <- l(dat$variable)
dat$upper1 <- u(dat$variable)
dat$variable1 <- v(dat$variable)

ggplot(dat,aes(x = variable)) + 
  geom_step(aes(y = variable1)) + 
  geom_ribbon(aes(ymin = upper1,ymax = lower1),alpha = 0.2)

可以让您从图中获取生成的数据，然后您可以过度绘制您喜欢的内容。

此图表显示：

红色=原装丝带

• 蓝色 = 采用原始 CI 向量并应用于 ecdf 曲线
• 绿色 = 计算上下序列和绘图的 ecdf

ggplot_build()

Answer 3

首先，我加载必要的包并生成虚拟数据。

    g<-ggplot(dat, aes(x = variable)) +
      geom_step(stat = "ecdf") +
      geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.5, fill="red") 

    inside<-ggplot_build(g)
    matched<-merge(inside$data[[1]],data.frame(x=dat$variable,dat$lower,dat$upper),by=("x"))

    g + 
      geom_ribbon(data=matched, aes(x = x, 
                                      ymin = y + dat.upper-x,
                                      ymax = y - x + dat.lower), 
                    alpha=0.5, fill="blue") +
      geom_ribbon(data=matched, aes(x = x, 
                                      ymin = ecdf(dat.lower)(x),
                                      ymax = ecdf(dat.upper)(x)), 
                    alpha=0.5, fill="green")

然后我生成一个图表：

library(reprex)
library(tidyverse)

set.seed(118)
data1 <- tibble(a = rnorm(100, mean = 20, sd = 3),
                ID = "a",
                gr = "A") %>% 
  bind_rows(tibble(a = rnorm(100, mean = 30, sd = 3),
                   ID = "z",
                   gr = "A")) %>% 
  bind_rows(tibble(a = rnorm(100, mean = 25, sd = 3),
                   ID = "a",
                   gr = "B")) %>% 
  bind_rows(tibble(a = rnorm(100, mean = 40, sd = 3),
                   ID = "z",
                   gr = "B")) %>% 
  mutate(CI_low = a - 3,
         CI_high = a + 3)

生成图表后，您可以预览其数据以及计算出的 ecdf：

data1 %>% 
  ggplot(aes(x = a, color = gr))+
  geom_step(stat = "ecdf")+
  facet_wrap(~ ID)

下一步是将分组变量的值分配给个人美学（在我的例子中是颜色）和变量组（方面的名称）。这必须手动完成，特别注意确保分配正确。

ggplot_build(last_plot())$data[[1]] %>% 
  as_tibble()
#> # A tibble: 408 × 10
#>    colour      y      x  ecdf flipped_aes PANEL group linewidth linetype alpha
#>    <chr>   <dbl>  <dbl> <dbl> <lgl>       <fct> <int>     <dbl>    <dbl> <lgl>
#>  1 #F8766D  0    -Inf    0    FALSE       1         1       0.5        1 NA   
#>  2 #F8766D  0.09   15.0  0.09 FALSE       1         1       0.5        1 NA   
#>  3 #F8766D  0.55   20.5  0.55 FALSE       1         1       0.5        1 NA   
#>  4 #F8766D  0.51   20.0  0.51 FALSE       1         1       0.5        1 NA   
#>  5 #F8766D  0.6    20.9  0.6  FALSE       1         1       0.5        1 NA   
#>  6 #F8766D  0.08   14.9  0.08 FALSE       1         1       0.5        1 NA   
#>  7 #F8766D  0.19   16.9  0.19 FALSE       1         1       0.5        1 NA   
#>  8 #F8766D  0.47   19.7  0.47 FALSE       1         1       0.5        1 NA   
#>  9 #F8766D  0.73   22.0  0.73 FALSE       1         1       0.5        1 NA   
#> 10 #F8766D  0.34   18.5  0.34 FALSE       1         1       0.5        1 NA   
#> # ℹ 398 more rows

下一步，我将数据与原始数据合并并生成一个新图表：

graph_data <- ggplot_build(last_plot())$data[[1]] %>% 
  as_tibble() %>% 
  mutate(gr = case_match(group,
                         1 ~ "A",
                         2 ~ "B")) %>% 
  mutate(ID = case_match(PANEL,
                         "1" ~ "a",
                         "2" ~ "z")) %>% 
  select(a = x, ecdf, gr, ID)

创建于 2024-07-25，使用

reprex v2.1.1