Можно ли воссоздать функциональность графика Bayesplot «mcmc_areas» в ggplot в R
Существует пакет, поддерживаемый Стэном, под названием bayesplot, который может создавать хорошие графики площади плотности с площадью под кривыми плотности, разделенными на основе интервалов достоверности для выборок апостериорных параметров, полученных через MCMC. В результате получается график, который выглядит следующим образом:
Я хочу создать аналогичный стиль графика с учетом одномерных списков выборочных значений с использованием ggplot, в который я могу передать любой общий список значений без необходимости подгонки Стэна и т. д. Кто-нибудь знает, как это сделать? Часть плотности понятна через geom_density, но я борюсь с разделением заполнения.
Ответы 2
Воспользуйтесь пакетом 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
Вот функция, которая генерирует график, похожий на 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")
Это отлично, спасибо, стоит отметить одну маленькую вещь: mcmc_areas использует интервалы с равными хвостами, а не максимальную плотность, на мгновение был сбит с толку тем, почему выходные данные немного отличаются, возможно, стоит изменить решение для использования eti (доступно в bayestestR package) или, по крайней мере, обратите внимание, что это так.
У меня остался один вопрос: как сделать так, чтобы гребни не перекрывались, это кажется вам нормальным, но когда я создаю графики с помощью вашей функции, гребни имеют тенденцию перекрывать друг друга и казаться намного больше?
Смотрите обновленный ответ. Функция теперь предоставляет выбор ETI или HDI и использует ETI по умолчанию. Я также показал, как использовать аргумент scale для управления перекрытием линий хребта.
Я получаю следующую ошибку с вашим примером: Joining, by = "name" Erreur : Problem with `filter()` input `..1`. ℹ Input `..1` is `between(x, CI_low, CI_high)`. x `left` must be length 1 ℹ The error occurred in group 1: name = "a".
Было бы полезно, если бы вы могли предоставить примеры данных для работы вместе с вашим вопросом. Это дает всем ответам общий набор данных для работы, что упрощает сравнение решений.