plot density plots with confidence intervals of 95% in R

Question

I'm trying draw multiple density plots in one plot for comparison porpuses. I wanted them to have their confidence interval of 95% like in the following figure. I'm working with ggplot2 and my df is a long df of observations for a certain location that I would like to compare for different time intervals.

I've done some experimentation following this example but I don't have the coding knowledge to achieve what I want. What i managed to do so far:

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]])

Help would be much apreaciated.

To add another density plot to that image. I tried to add another plot but i get an error message saying iIcan't do that to ggplto() object. And also a mean line. — noob_researcher, Oct 27 '21 at 17:21
Several posts that should have you covered: https://stackoverflow.com/q/4542438/5325862, https://stackoverflow.com/q/41971150/5325862, https://stackoverflow.com/q/64227409/5325862, https://stackoverflow.com/q/49807993/5325862 — camille, Oct 27 '21 at 17:49

score 0 · Answer 1 · answered Oct 27 '21 at 20:05

This is obviously on a different set of data, but this is roughly that plot with data from 2 t distributions. I've included the data generation in case it is of use.

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))

plot density plots with confidence intervals of 95% in R

1 Answers1