color_under_density_facets.R

library(dplyr)
library(tidyr)
library(magrittr)
library(bayesplot)

plot_areas <- function(d, n, mu, s, cutoff, left = TRUE, h, v) {
  
  if (left) {
    Pr <- apply(d, 2, function (x) mean(x < cutoff))
  } else {
    Pr <- apply(d, 2, function (x) mean(x >= cutoff))
  }
  key <- colnames(d)
  d <- gather(d)
  d %<>% 
    group_by(key) %>%
    mutate(x = density(value, bw = "SJ", n = n)$x, 
           y = density(value, bw = "SJ", n = n)$y,
           x_pos = mean(x),
           y_pos = mean(y))
  if (left) {
    d <- mutate(d, area = x < cutoff)
  } else {
    d <- mutate(d, area = x >= cutoff)
  }
  
  label <- paste("Pr =", as.character(round(Pr, 2)))
  dat_text <- data.frame(label = label, key = key)
  
  p <- ggplot(d, aes(x, y)) +
    geom_ribbon(aes(ymin = 0, ymax = y, fill = area), alpha = 1/2) +
    facet_wrap(~key) +
    geom_vline(xintercept = cutoff, size = 0.2, linetype = 2) +
    ylab("") +
    geom_text(
      data    = dat_text,
      mapping = aes(x = -Inf, y = -Inf, label = label),
      hjust   = h,
      vjust   = v
    ) + 
    bayesplot::theme_default() + 
    theme(legend.position="none",
          axis.title.y = element_blank(),
          axis.text.y  = element_blank(),
          axis.ticks.y = element_blank(),
          axis.line.y  = element_blank()) 
  return(p)
}

n <- 1e3
mu <- c(0.69, 0.76, 0.81, 0.78)
s <- c(0.1, 0.12, 0.15, 0.17)
cutoff <- 0.6
d <- tibble(A = rnorm(n, mu[1], s[1]),
            B = rnorm(n, mu[2], s[2]),
            C = rnorm(n, mu[3], s[3]),
            D = rnorm(n, mu[4], s[4]))
p <- plot_areas(d, n, mu, s, cutoff, TRUE, -0.5, -3) + xlab("Hazard Ratio")
p