在 R 中以 95% 的置信区间绘制密度图

Posted 2023-02-16

【中文标题】在 R 中以 95% 的置信区间绘制密度图

【英文标题】：plot density plots with confidence intervals of 95% in R

【发布时间】：2021-12-12 23:33:45

【问题描述】：

我正在尝试在一个图中绘制多个密度图以进行比较。我希望他们的置信区间为 95%，如下图所示。我正在使用ggplot2，我的 df 是某个位置的长 df 观察值，我想比较不同的时间间隔。

我在example 之后做了一些实验，但我没有实现我想要的编码知识。 到目前为止，我设法做到了：

library(magrittr)
library(ggplot2)
library(dplyr)

build_object <- ggplot_build(
  ggplot(data=ex_long, aes(x=val)) + geom_density())

plot_credible_interval <- function(
  gg_density,  # ggplot object that has geom_density
  bound_left,
  bound_right
) 
  build_object <- ggplot_build(gg_density)
  x_dens <- build_object$data[[1]]$x
  y_dens <- build_object$data[[1]]$y
  
  index_left <- min(which(x_dens >= bound_left))
  index_right <- max(which(x_dens <= bound_right))
  
  gg_density + geom_area(
    data=data.frame(
      x=x_dens[index_left:index_right],
      y=y_dens[index_left:index_right]), 
    aes(x=x,y=y),
    fill="grey",
    alpha=0.6)


gg_density <- ggplot(data=ex_long, aes(x=val)) + 
  geom_density()
gg_density %>% plot_credible_interval(tab$q2.5[[40]], tab$q97.5[[40]])

非常感谢您的帮助。

【问题讨论】：

那不是你想要的呢？

为该图像添加另一个密度图。我试图添加另一个情节，但我收到一条错误消息，说 iIcan't do that to ggplto() object。还有一条平均线。

应该让您了解的几个帖子：***.com/q/4542438/5325862、***.com/q/41971150/5325862、***.com/q/64227409/5325862、***.com/q/49807993/5325862

【参考方案1】：

这显然是在一组不同的数据上，但这大致是来自 2 个 t 分布的数据的图。我已经包含了数据生成以备不时之需。

library(tidyverse)

x1 <- seq(-5, 5, by = 0.1)
t_dist1 <- data.frame(x = x1,
                     y = dt(x1, df = 3),
                     dist = "dist1")
x2 <- seq(-5, 5, by = 0.1)
t_dist2 <- data.frame(x = x2,
                      y = dt(x2, df = 3),
                      dist = "dist2")

t_data = rbind(t_dist1, t_dist2) %>%
  mutate(x = case_when(
    dist == "dist2" ~ x + 1,
    TRUE ~ x
  ))

p <- ggplot(data = t_data,
            aes(x = x,
                y = y )) +
  geom_line(aes(color = dist))

plot_data <- as.data.frame(ggplot_build(p)$data)

bottom <- data.frame(plot_data) %>%
  mutate(dist = case_when(
    group == 1 ~ "dist1",
    group == 2 ~ "dist2"
  )) %>%
  group_by(dist) %>%
  slice_head(n = ceiling(nrow(.) * 0.1)) %>% 
  ungroup()

top <- data.frame(plot_data) %>%
  mutate(dist = case_when(
    group == 1 ~ "dist1",
    group == 2 ~ "dist2"
  )) %>%
  group_by(dist) %>%
  slice_tail(n = ceiling(nrow(.) * 0.1)) %>%
  ungroup()

segments <- t_data %>%
  group_by(dist) %>%
  summarise(x = mean(x),
            y = max(y))

p + geom_area(data = bottom,
              aes(x = x,
                  y = y,
                  fill = dist),
              alpha = 0.25,
              position = "identity") +
  geom_area(data = top,
            aes(x = x,
                y = y,
                fill = dist),
            alpha = 0.25,
            position = "identity") +
  geom_segment(data = segments,
               aes(x = x,
                   y = 0,
                   xend = x,
                   yend = y,
                   color = dist,
                   linetype = dist)) +
  scale_color_manual(values = c("red", "blue")) +
  scale_linetype_manual(values = c("dashed", "dashed"),
                        labels = NULL) +
  ylab("Density") +
  xlab("\U03B2 for AQIv") +
  guides(color = guide_legend(title = "p.d.f \U03B2",
                              title.position = "right",
                              labels = NULL),
         linetype = guide_legend(title = "Mean \U03B2",
                                 title.position = "right",
                                 labels = NULL,
                                 override.aes = list(color = c("red", "blue"))),
         fill = guide_legend(title = "Rej. area \U03B1 = 0.05",
                             title.position = "right",
                             labels = NULL)) +
  annotate(geom = "text",
           x = c(-4.75, -4),
           y = 0.35,
           label = c("RK", "OK")) +
  theme(panel.background = element_blank(),
        panel.border = element_rect(fill = NA,
                                    color = "black"),
        legend.position = c(0.2, 0.7),
        legend.key = element_blank(),
        legend.direction = "horizontal",
        legend.text = element_blank(),
        legend.title = element_text(size = 8))