是否有可能在 R 中的 ggplot 中重新创建 bayesplot 的 "mcmc_areas" 图的功能

Question

Stan 支持一个名为 bayesplot 的包，它可以生成漂亮的密度区域图，其中密度曲线下的区域根据通过 MCMC 绘制的后验参数样本的可信区间进行分区，这导致情节如下所示：

我希望使用 ggplot 在给定一维采样值列表的情况下制作类似样式的绘图，我可以将任何通用值列表传递给它而不必是 Stan 拟合等。有谁知道该怎么做这个？密度部分通过 geom_density 很清楚，但我正在为填充分区而苦苦挣扎。

Answer 1

使用ggridges包：

library(tidyverse)
library(ggridges)

tibble(data_1, data_2, data_3) %>%
  pivot_longer(everything()) %>%
  ggplot(aes(x = value, y = name, group = name)) +
  geom_density_ridges()

数据：

set.seed(123)
n <- 15
data_1 <- rnorm(n)
data_2 <- data_1 - 1
data_3 <- data_1 + 2

Answer 2

这是一个生成类似于 bayesplot::mcmc_areas 的绘图的函数。它绘制可信区间（默认等尾，或最高密度），可选择设置区间的概率宽度：

library(tidyverse)
library(ggridges)
library(bayestestR)
theme_set(theme_classic(base_size=15))

# Create ridgeplots with credible intervals
# ARGUMENTS
# data A data frame
# FUN A function that calculates credible intervals
# ci The width of the credible interval
# ... For passing optional arguments to geom_ridgeline. 
#      For example, change the scale parameter to control overlap of ridge lines. 
#      geom_ridgeline's default is scale=1.
plot_density_ridge = function(data, FUN=c("eti", "hdi"), ci=0.89, ...) {
  
  # Determine whether to use eti or hdi function
  FUN = match.arg(FUN)
  FUN = match.fun(FUN)
  
  # Get kernel density estimate as a data frame
  dens = map_df(data, ~ {
    d = density(.x, na.rm=TRUE)
    tibble(x=d$x, y=d$y)
  }, .id="name")
  
  # Set relative width of median line
  e = diff(range(dens$x)) * 0.006
  
  # Get credible interval width and median
  cred.int = data %>% 
    pivot_longer(cols=everything()) %>% 
    group_by(name) %>% 
    summarise(CI=list(FUN(value, ci=ci)),
              m=median(value, na.rm=TRUE)) %>% 
    unnest_wider(CI)
  
  dens %>% 
    left_join(cred.int) %>% 
    ggplot(aes(y=name, x=x, height=y)) +
      geom_vline(xintercept=0, colour="grey70") +
      geom_ridgeline(data= . %>% group_by(name) %>%
                       filter(between(x, CI_low, CI_high)),
                     fill=hcl(230,25,85), ...) +
      geom_ridgeline(data=. %>% group_by(name) %>% 
                       filter(between(x, m - e, m + e)),
                     fill=hcl(240,30,60), ...) +
      geom_ridgeline(fill=NA, ...) + 
      geom_ridgeline(fill=NA, aes(height=0), ...) +
      labs(y=NULL, x=NULL)
  
}

现在我们来试试这个功能

# Fake data
set.seed(2)
d = data.frame(a = rnorm(1000, 0.6, 1),
               b = rnorm(1000, 1.3, 0.5),
               c = rnorm(1000, -1.2, 0.7))

plot_density_ridge(d)
plot_density_ridge(d, ci=0.5, scale=1.5)
plot_density_ridge(iris %>% select(-Species)) 
plot_density_ridge(iris %>% select(-Species), FUN="hdi")

是否有可能在 R 中的 ggplot 中重新创建 bayesplot 的 "mcmc_areas" 图的功能

Is it possible to recreate the functionality of bayesplot's "mcmc_areas" plot in ggplot in R

r

bayesian

ggplot2

stan