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

Question

在阅读了不同的帖子后，我发现了如何向密度图添加一条均值的 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？

Answer 1

使用来自链接的数据，您可以像这样计算均值、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_segment 结合起来。 s:

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)