同一图中不同密度曲线（分组因子）下的阴影面积

Question

我正在尝试绘制一条密度线，并且只想对与 x 轴的 95% 相关联的区域进行着色或填充。我正在尝试遵循所附答案中给出的答案，但是当我们使用分组因子同时绘制多个分布时，他们都没有谈论阴影区域。在这种情况下，分组因子是不同的中心电极（“Fz”、“Cz”、“Pz”）。我试图将类似于最高密度间隔的东西可视化，或者曲线下面积介于 5 和 95 之间。

Area Under Curve AUC by Group

我的数据如下所示：

> head(dframe1)
         x            y Electrode
1 1.571296 0.0001474116        Fz
2 1.576496 0.0001487649        Fz
3 1.581697 0.0001497564        Fz
4 1.586897 0.0001504074        Fz
5 1.592098 0.0001507446        Fz
6 1.597298 0.0001507776        Fz

目前我用来在 ggplot 中按组绘制分布的代码如下所示：

p1 <- ggplot(data = dframe1, mapping = aes(x = x, y = y)) +
  geom_density_line(stat = "identity", size=.5, alpha=0.3, aes(color=Electrode, fill=Electrode)) +
  scale_fill_discrete(breaks=c("Fz","Cz","Pz")) +
  guides(colour = FALSE) +
  geom_vline(xintercept = 0) +
  xlab("values") +
  xlim(-2, 10) +
  ylab("density") +
  ylim(0, .7) +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=16),
        plot.title = element_text(size=18)) +
  labs(title="Interval")

我画了一些类似于我正在寻找的东西：

当然我可以使用bayestestR HDI标准输出，但我更喜欢ggplot的美感和灵活性。

任何帮助将不胜感激。

Answer 1

由于您需要对密度曲线进行一些数学运算来计算出 95% 区间的位置，因此最好在 在 ggplot 之外执行此操作。我经常发现人们遇到问题是因为他们试图让 ggplot 做太多的数据争论和总结。通常更容易计算出您想要绘制的内容，然后绘制它。

在您的情况下，您的 x 和 y 坐标已经代表密度。对于每个电极，您只需要创建一个逻辑向量，告诉您密度的积分何时介于 0.025 和 0.975 之间，这样您就可以轻松地从 95% 置信区间中提取子集。您可以使用 split-aplly-bind 方法来做到这一点，如下所示：

densdf <- do.call(rbind, lapply(split(dframe1, dframe1$Electrode), function(z)
{
  integ <- cumsum(z$y * mean(diff(z$x)))
  CI <-  integ > 0.025 & integ < 0.975
  data.frame(x = z$x, y = z$y, Electrode = z$Electrode[1], CI = CI)
}))

现在我们准备好绘图了：

ggplot(data = densdf, mapping = aes(x = x, y = y)) +
  geom_area(data = densdf[densdf$CI,], 
            aes(fill = Electrode, color = Electrode),
            outline.type = "full", alpha = 0.3, size = 1) +
  geom_line(aes(color = Electrode), size = 1) +
  scale_fill_discrete(breaks = c("Fz", "Cz", "Pz")) +
  guides(colour = FALSE) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  lims(x = c(-2, 10), y = c(0, 0.7)) +
  labs(title = "Interval", x = "values", y = "density") +
  theme_bw() +
  theme(axis.text  = element_text(size = 12),
        axis.title = element_text(size = 16),
        plot.title = element_text(size = 18))