目标是展示各个ROI后验分布，95%HDI以及Mode值。在Kruschke的书中，他就提供了可以直接生成下图的函数。

当有多个ROI的后验分布需要展示时，山脊图就成了不错的选择。ggridges的官方教程(https://wilkelab.org/ggridges/articles/introduction.html)一看就会不多做介绍，这里只记录教程里没怎么提及的五毛特效。

1. 添加线段显示分布的统计值(Mode)

2. 展示分布的95%HDI区间

3. 根据统计值进行排序

4. 给山脊图添加阴影

最后效果如下：

配色使用了Kruschke的配色: skyblue

1.如何计算95%HDI

平时常用的95%的ETI(equal-tailed interval)在其界限的两侧有2.5%的分布,分别代表了2.5百分位数和97.5百分位数。在对称分布中，ETI和HDI(higest density interval)是一样的，但在非对称分布中则不然，例如下图的伽马分布，标出了95%的HDI和95%的ETI，显然使用95%HDI比95%ETI更准确地反映了分布的可信度。

95%HDI可以用Kruschke给的代码计算，也可以使用HDInterval中的hdi函数计算，后者支持的object类型很多，可以直接使用rjags得到的mcmc.list作为输入，默认计算的是95%HDI。

2.如何计算mcmc后验分布的mode值

在最大后验估计（MAP）中，我们使用后验分布的最大值来推导出我们感兴趣的点估计。一个分布的最大值也被称为它的 "mode"（假设是unimodal分布)。用mean，median还是mode作为后验分布的点估计可能依据情况而定。在Kruschke的书中他推荐了mode，但也说median相对mode会更稳定，在他的博客中他直接推荐的是mode。也有观点认为在某些条件下mean也是可行的。

将Kruschke提供的代码做出函数，可计算一个vector的MAP。

get_mode <- function(paramSampleVec){mcmcDensity = density(paramSampleVec)mo = mcmcDensity$x[which.max(mcmcDensity$y)]return(mo)
}

3. 数据准备

数据需要做出key+value的形式，使用gather转换为长表格形式(如下)，之后就可以用来作图了。

df.para=as.data.frame(all_post.para) %>% gather(key='Region',value = 'Value')

4.生成一个基础版的山脊图

ggplot(df.para)+geom_density_ridges(aes(x = Value, y = Region)

5.根据教程内容添加几条线

很简单，就是加上quantile_lines，这里显示0.025和0.975的quantiles，也就是之前提到的95%ETI。

ggplot(df.para)+stat_density_ridges(aes(x = Value, y = Region),quantile_lines = TRUE, quantiles = c(0.025, 0.975))

6.把添加的线条换成我们想要的mode和95%HDI

如果需要显示自定义的统计值，比如mode和95%HDI的话，需添加一个quantile_fun，该函数返回相应统计值即可。结合计算mode和hdi的方法，建立一个函数，同时返回95% HDI lower，mode和HID upper三个值。

mode_hdi <- function(x,...){hdi_bound=hdi(x)mode=get_mode(x)return(c(hdi_bound['lower'],mode,hdi_bound['upper']))
}

再将这个函数放到ggridges中即可，可用vline_linetype设置线条类型。

ggplot(df.para)+stat_density_ridges(aes(x = Value, y = Region),quantile_lines = TRUE, quantile_fun = mode_hdi,vline_linetype = 2)

7. 给95%HDI加上颜色

在aes中加入fill = factor(stat(quantile))

之后将stat_density_ridges中加入geom = "density_ridges_gradient"

不设置geom，也可以直接使用geom_density_ridges_gradient函数。

此时fill的内容就被quantile_fun的三条线分为了四个部分，用scale_fill_manual添加四个自选的颜色，中间两个便是95%HDI。

ggplot(df.para)+stat_density_ridges(aes(x = Value, y = Region,fill = factor(stat(quantile))),geom = "density_ridges_gradient",quantile_lines = TRUE, quantile_fun = mode_hdi,vline_linetype = 2)+scale_fill_manual(values=c("skyblue4", "skyblue", "skyblue","skyblue4"), guide = "none")

8. 根据mode排序

在数据准备时，用fct_reorder将Region根据mode值排序即可，其中的.fun使用的是计算mode的函数。依此类推，设置.fun函数，就可以根据自定义的统计量进行排序，例如95%HDI的宽度，mean，median值等。

df.para=as.data.frame(all_post.para) %>% gather(key='Region',value = 'Value') %>% mutate(Region=fct_reorder(Region,Value,.fun=get_mode,.desc = TRUE))

用该数据作图如下，根据mode排序完成：

9. 替换ROI的缩写

可以准备一个缩写的对照表如下：

在数据准备时，用rename_all和str_replace_all根据对照表将缩写批量替换。

df.para=as.data.frame(all_post.para) %>% rename_all(~stringr::str_replace(.,"L_|R_","")) %>% # Rename labelsrename_all(~stringr::str_replace_all(.,setNames(full_label,paste0(orig_label,'$')))) %>% # add $ to get exact match gather(key='Region',value = 'Value') %>% mutate(Region=fct_reorder(Region,Value,.fun=get_mode,.desc = TRUE))

用该数据作图如下，替换label完成：

额外内容

1. 添加阴影和垂直线

在山脊图的代码加入geom_rect添加一个有颜色的矩形即可，添加垂直线也是一样的原理。

geom_rect(data = data.frame(x = 1),xmin = -0.2, xmax = -0.05, ymin = -Inf, ymax = Inf,alpha = 0.4, fill = "gray")+
geom_vline(xintercept = -0.1)

2. 可以考虑将mode用ggseg画出来，并跟山脊图拼在一起。

->ggseg教程<-

3. bayestestR和bayesplot等工具包都有一些展示后验分布的函数可供使用。

bayestestR

bayesplot