将 vline 添加到 geom_density 和平均 R 的阴影置信区间

Posted 2023-02-16

【中文标题】将 vline 添加到 geom_density 和平均 R 的阴影置信区间

【英文标题】：Add vline to geom_density and shade confidence interval of mean R

【发布时间】：2017-06-17 16:34:33

【问题描述】：

在阅读了不同的帖子后，我发现了如何在密度图中添加一条均值 vline，如 here 所示。 使用上述链接中提供的数据：

1) 如何使用 geom_ribbon 在平均值周围添加 95% 的置信区间？ CI 可以计算为

#computation of the standard error of the mean
sem<-sd(x)/sqrt(length(x))
#95% confidence intervals of the mean
c(mean(x)-2*sem,mean(x)+2*sem)

2) 如何将 vline 限制在曲线下的区域？您将在下图中看到曲线外的 vline 绘图。

可以在https://www.dropbox.com/s/bvvfdpgekbjyjh0/test.csv?dl=0找到与我的实际问题非常接近的示例数据

更新

使用上面链接中的真实数据，我使用@beetroot的答案尝试了以下内容。

# Find the mean of each group
dat=me
library(dplyr)
library(plyr)
cdat <- ddply(data,.(direction,cond), summarise, rating.mean=mean(rating,na.rm=T))# summarize by season and variable
cdat

#ggplot
p=ggplot(data,aes(x = rating)) + 
  geom_density(aes(colour = cond),size=1.3,adjust=4)+
  facet_grid(.~direction, scales="free")+
  xlab(NULL) + ylab("Density")
p=p+coord_cartesian(xlim = c(0, 130))+scale_color_manual(name="",values=c("blue","#00BA38","#F8766D"))+
  scale_fill_manual(values=c("blue", "#00BA38", "#F8766D"))+
  theme(legend.title = element_text(colour="black", size=15, face="plain"))+
  theme(legend.text = element_text(colour="black", size = 15, face = "plain"))+
  theme(title = red.bold.italic.text, axis.title = red.bold.italic.text)+
  theme(strip.text.x = element_text(size=20, color="black",face="plain"))+ # facet labels
  ggtitle("SAMPLE A") +theme(plot.title = element_text(size = 20, face = "bold"))+
    theme(axis.text = blue.bold.italic.16.text)+ theme(legend.position = "none")+
  geom_vline(data=cdat, aes(xintercept=rating.mean, color=cond),linetype="dotted",size=1)
p

## implementing @beetroot's code to restrict lines under the curve and shade CIs around the mean
# I will use ddply for mean and CIs
cdat <- ddply(data,.(direction,cond), summarise, rating.mean=mean(rating,na.rm=T),
              sem = sd(rating,na.rm=T)/sqrt(length(rating)),
              ci.low = mean(rating,na.rm=T) - 2*sem,
              ci.upp = mean(rating,na.rm=T) + 2*sem)# summarize by direction and variable


#In order to limit the lines to the outline of the curves you first need to find out which y values
#of the curves correspond to the means, e.g. by accessing the density values with ggplot_build and 
#using approx:

   cdat.dens <- ggplot_build(ggplot(data, aes(x=rating, colour=cond)) +
                              facet_grid(.~direction, scales="free")+
                              geom_density(aes(colour = cond),size=1.3,adjust=4))$data[[1]] %>%
  mutate(cond = ifelse(group==1, "A",
                       ifelse(group==2, "B","C"))) %>%
  left_join(cdat) %>%
  select(y, x, cond, rating.mean, sem, ci.low, ci.upp) %>%
  group_by(cond) %>%
  mutate(dens.mean = approx(x, y, xout = rating.mean)[[2]],
         dens.cilow = approx(x, y, xout = ci.low)[[2]],
         dens.ciupp = approx(x, y, xout = ci.upp)[[2]]) %>%
  select(-y, -x) %>%
  slice(1)

 cdat.dens

#---
 #You can then combine everything with various geom_segments:

   ggplot(data, aes(x=rating, colour=cond)) +
   geom_density(data = data, aes(x = rating, colour = cond),size=1.3,adjust=4) +facet_grid(.~direction, scales="free")+
   geom_segment(data = cdat.dens, aes(x = rating.mean, xend = rating.mean, y = 0, yend = dens.mean, colour = cond),
                linetype = "dashed", size = 1) +
   geom_segment(data = cdat.dens, aes(x = ci.low, xend = ci.low, y = 0, yend = dens.cilow, colour = cond),
                linetype = "dotted", size = 1) +
   geom_segment(data = cdat.dens, aes(x = ci.upp, xend = ci.upp, y = 0, yend = dens.ciupp, colour = cond),
                linetype = "dotted", size = 1)

给出这个：

您会注意到平均值和 CI 没有像原始图中那样对齐。 @beetroot 我做错了什么？

【问题讨论】：

我会使用 geom_rect。将 ymin 和 ymax 设置为 +/-Inf

为了限制曲线，您必须预先计算密度。

另见***.com/questions/12429333/…

我刚刚意识到我没有跟进您的问题..您最终弄清楚了吗？

@beetroot 下面接受的解决方案对我的数据非常有效。感谢您跟进此事。

【参考方案1】：

如果您想绘制平均线而不构建绘图对象并且在绘图之前不操纵数据，您可以使用stat_summary()：

(
    ggplot(data = dat, aes(x = rating, colour = cond))
    + geom_density()
    + stat_summary(
        aes(y = rating, x = 0),
        geom = 'rect',
        fun.data = density_mean_line(dat$rating),
        key_glyph = "vline",
        size = 1
    )
)

给予：

地点：

density_mean_line = function(values) 
    values_range = range(values, na.rm=TRUE)
    function(x) 
        density_data = StatDensity$compute_group(
            data.frame(x=x),
            scales=list(
                x=scale_x_continuous(limits = values_range)
            )
        )
        mean_x = mean(x)
        data.frame(
            xmin=mean_x,
            xmax=mean_x,
            ymin=0,
            ymax=approx(density_data$x, density_data$density, xout=mean_x)$y
        )
    

而dat 被定义为erc 的答案：

set.seed(1234)
dat <- data.frame(cond = factor(rep(c("A","B"), each=200)), 
                  rating = c(rnorm(200),rnorm(200, mean=.8)))

此技术也可用于生成实心区域（与密度轮廓颜色相同）：

(
    ggplot(data = dat, aes(x = rating, colour = cond, group = cond))
    + stat_summary(
        aes(y = rating, x = 0, fill = cond),
        geom = 'rect',
        fun.data = density_ci(dat$rating),
        size=1
    )
    + stat_summary(
        aes(y = rating, x = 0),
        geom = 'rect',
        fun.data = density_mean_line(dat$rating),
        key_glyph = "vline",
        size = 0.5,
        color='grey20'
    )
    + geom_density()
)

地点：

density_ci = function(values, resolution=100) 
    values_range = range(values, na.rm=TRUE)
    function(x) 
        density_data = StatDensity$compute_group(
            data.frame(x=x),
            scales=list(
                x=scale_x_continuous(limits = values_range)
            )
        )
        mean_x = mean(x)
        sem = sd(x) / sqrt(length(x))
        ci_lower = mean_x - 1.96 * sem
        ci_upper = mean_x + 1.96 * sem

        x_values = seq(ci_lower, ci_upper, length.out=resolution)

        data.frame(
            xmin=x_values,
            xmax=x_values,
            ymin=rep(0, resolution),
            ymax=approx(density_data$x, density_data$density, xout=x_values)$y
        )
    

【参考方案2】：

使用链接中的数据，您可以像这样计算平均值，se和ci（我建议使用dplyr，plyr的继任者）：

set.seed(1234)
dat <- data.frame(cond = factor(rep(c("A","B"), each=200)), 
                  rating = c(rnorm(200),rnorm(200, mean=.8)))

library(ggplot2)
library(dplyr)
cdat <- dat %>%
  group_by(cond) %>%
  summarise(rating.mean = mean(rating),
            sem = sd(rating)/sqrt(length(rating)),
            ci.low = mean(rating) - 2*sem,
            ci.upp = mean(rating) + 2*sem)

为了将线条限制在曲线的轮廓上，您首先需要找出曲线的哪些 y 值对应于均值，例如通过ggplot_build 和approx 访问密度值：

cdat.dens <- ggplot_build(ggplot(dat, aes(x=rating, colour=cond)) + geom_density())$data[[1]] %>%
  mutate(cond = ifelse(group == 1, "A", "B")) %>%
  left_join(cdat) %>%
  select(y, x, cond, rating.mean, sem, ci.low, ci.upp) %>%
  group_by(cond) %>%
  mutate(dens.mean = approx(x, y, xout = rating.mean)[[2]],
         dens.cilow = approx(x, y, xout = ci.low)[[2]],
         dens.ciupp = approx(x, y, xout = ci.upp)[[2]]) %>%
  select(-y, -x) %>%
  slice(1)

> cdat.dens
Source: local data frame [2 x 8]
Groups: cond [2]

   cond rating.mean        sem     ci.low     ci.upp dens.mean dens.cilow dens.ciupp
  <chr>       <dbl>      <dbl>      <dbl>      <dbl>     <dbl>      <dbl>      <dbl>
1     A -0.05775928 0.07217200 -0.2021033 0.08658471 0.3865929   0.403623  0.3643583
2     B  0.87324927 0.07120697  0.7308353 1.01566320 0.3979347   0.381683  0.4096153

然后您可以将所有内容与各种geom_segments 结合起来：

ggplot() +
  geom_density(data = dat, aes(x = rating, colour = cond)) +
  geom_segment(data = cdat.dens, aes(x = rating.mean, xend = rating.mean, y = 0, yend = dens.mean, colour = cond),
             linetype = "dashed", size = 1) +
  geom_segment(data = cdat.dens, aes(x = ci.low, xend = ci.low, y = 0, yend = dens.cilow, colour = cond),
             linetype = "dotted", size = 1) +
  geom_segment(data = cdat.dens, aes(x = ci.upp, xend = ci.upp, y = 0, yend = dens.ciupp, colour = cond),
               linetype = "dotted", size = 1)

正如 Axeman 指出的那样，您可以根据 this answer 中的说明根据功能区区域创建多边形。

因此，对于您的数据，您可以像这样子集并添加其他行：

ribbon <- ggplot_build(ggplot(dat, aes(x=rating, colour=cond)) + geom_density())$data[[1]] %>%
  mutate(cond = ifelse(group == 1, "A", "B")) %>%
  left_join(cdat.dens) %>%
  group_by(cond) %>%
  filter(x >= ci.low & x <= ci.upp) %>%
  select(cond, x, y)

ribbon <- rbind(data.frame(cond = c("A", "B"), x = c(-0.2021033, 0.7308353), y = c(0, 0)), 
                as.data.frame(ribbon), 
                data.frame(cond = c("A", "B"), x = c(0.08658471, 1.01566320), y = c(0, 0)))

并将geom_polygon 添加到情节中：

ggplot() +
  geom_polygon(data = ribbon, aes(x = x, y = y, fill = cond), alpha = .5) +
  geom_density(data = dat, aes(x = rating, colour = cond)) +
  geom_segment(data = cdat.dens, aes(x = rating.mean, xend = rating.mean, y = 0, yend = dens.mean, colour = cond),
             linetype = "dashed", size = 1) +
  geom_segment(data = cdat.dens, aes(x = ci.low, xend = ci.low, y = 0, yend = dens.cilow, colour = cond),
             linetype = "dotted", size = 1) +
  geom_segment(data = cdat.dens, aes(x = ci.upp, xend = ci.upp, y = 0, yend = dens.ciupp, colour = cond),
               linetype = "dotted", size = 1)

---

这是为您的真实数据改编的代码。合并两个组而不是一个组有点棘手：

cdat <- dat %>%
  group_by(direction, cond) %>%
  summarise(rating.mean = mean(rating, na.rm = TRUE),
            sem = sd(rating, na.rm = TRUE)/sqrt(length(rating)),
            ci.low = mean(rating, na.rm = TRUE) - 2*sem,
            ci.upp = mean(rating, na.rm = TRUE) + 2*sem)

cdat.dens <- ggplot_build(ggplot(dat, aes(x=rating, colour=interaction(direction, cond))) + geom_density())$data[[1]] %>%
  mutate(cond = ifelse((group == 1 | group == 2 | group == 3 | group == 4), "A",
                        ifelse((group == 5 | group == 6 | group == 7 | group == 8), "B", "C")),
         direction = ifelse((group == 1 | group == 5 | group == 9), "EAST",
                            ifelse((group == 2 | group == 6 | group == 10), "NORTH",
                                   ifelse((group == 3 | group == 7 | group == 11), "SOUTH", "WEST")))) %>%
  left_join(cdat) %>%
  select(y, x, cond, direction, rating.mean, sem, ci.low, ci.upp) %>%
  group_by(cond, direction) %>%
  mutate(dens.mean = approx(x, y, xout = rating.mean)[[2]],
         dens.cilow = approx(x, y, xout = ci.low)[[2]],
         dens.ciupp = approx(x, y, xout = ci.upp)[[2]]) %>%
  select(-y, -x) %>%
  slice(1)

ggplot() +
  geom_density(data = dat, aes(x = rating, colour = cond)) +
  geom_segment(data = cdat.dens, aes(x = rating.mean, xend = rating.mean, y = 0, yend = dens.mean, colour = cond),
               linetype = "dashed", size = 1) +
  geom_segment(data = cdat.dens, aes(x = ci.low, xend = ci.low, y = 0, yend = dens.cilow, colour = cond),
               linetype = "dotted", size = 1) +
  geom_segment(data = cdat.dens, aes(x = ci.upp, xend = ci.upp, y = 0, yend = dens.ciupp, colour = cond),
               linetype = "dotted", size = 1) +
  facet_wrap(~direction)

【讨论】：

哇！这看起来很神奇。我会在我的数据上尝试一下，然后告诉你。再次感谢。

请问您可以在我上面提供的数据上使用您的代码吗？我的实际数据有condA,B,C 和一个附加变量direction，即east,west,south,north。最终目标是使用facet_grid(.~direction, scales="free") 显示密度，我绝对可以做到。我还有更多的小组要进行分析，但是您对上面链接中新提供的数据test.csv 的回答应该让我开始。谢谢。

请在您的问题中包含您的数据的dput()，并说明您所做的工作，并说明您在应用我的回答中的代码时遇到问题的地方。

我刚刚添加了一个修改后的代码，显示了到目前为止我用示例图与初始图相比所做的工作。数据太大，无法生成dput，但可以在上面的 libk 中找到。感谢您为分析提供更多信息。

当我运行cdat <- dat %>% group_by(direction, cond) %>% summarise(rating.mean = mean(rating, na.rm = TRUE), sem = sd(rating, na.rm = TRUE)/sqrt(length(rating)), ci.low = mean(rating, na.rm = TRUE) - 2*sem, ci.upp = mean(rating, na.rm = TRUE) + 2*sem) 时，我得到的只是1 Obs of 4 variables。然后当我运行cdat.dens 我得到Error: No common variables. Please specify by param. 是否汇总忽略direction or cond？