简介
最近科研过程中，需要可视化系统剩余使用寿命（Remaining Useful Lifetime, RUL）的分布。刚开始的想法是使用这篇推文：复现 sci 顶刊中的 3D 密度函数图。
后来发现在部署 R 包中会存在一个 BUG，并且图片的大小的美观度不算最佳。
所以这次尝试使用峰峦图来可视化结果。2021年的时候写了一篇 ggridges 包的教程：ggridges包—峰峦图详细介绍。写的比较详细，感兴趣的读者可以翻阅。
加载相关 R 包：
library(tidyverse) library(ggridges) library(viridis)
数据介绍
这里用一个模拟的数据作为例子。实际过程中，在经过贝叶斯后验抽样也会得到类似这样的结果。
假设绘制未来 5 天（y <- c(1,6,11,16,21)）的 RUL 分布，以正态分布作为模拟分布 N(mean1,std1)。
mean1 = seq(21,1,length.out =5) # 刻画不同时间对应密度函数的均值 std1 = seq(5,3,length.out =5) # 刻画不同时间对应密度函数的标准差 len = 10000 #分布中离散数据的个数 y <- c(1,6,11,16,21) # 未来时间 z = matrix(NA,len,5) # 存储抽样得到的后验 for(i in 1:length(mean1)){ z[,i] = rnorm(len,mean1[i],std1[i]) } # 合并数据 dat = as.data.frame(cbind(y,z)) colnames(dat) = c("Cur_time",paste("z",1:5,sep = "")) head(dat)
该数据中 z1:z5 表示未来 5 天的后验结果。这里的数据不算是整洁数据。需要通过 pivot_longer() 进行数据整理。
dat %>% pivot_longer(`z1`:`z5`, names_to = "Class", values_to = "Value" ) -> new_dat
我们将 z1:z5 名称存储到 "Class" 中，数值存储到 "Value" 中。转化后的数据如下所示：
初级版本
根据以上数据 new_dat，我们先使用 ggridges[1] 包中的 geom_density_ridges() 进行简单的可视化。
ggplot(new_dat, aes(x = Value, y = Class, fill = Class)) + geom_density_ridges()
基本就是小编想要的结果，x 轴表示 RUL（数值变量），y 轴表示未来某天（分类变量）。之后需要在这个基础上按照自身需求做一些调整，可参考：ggridges包—峰峦图详细介绍。
接下来，小编给出几种自己常用的方式。
进阶版本
细节处理
小编通常使用 theme_bw() 主题，并配合 theme(panel.grid = element_blank())使用。填充颜色经常使用 viridis 包中的 scale_fill_viridis()。
关于 y 轴标签名字的修改可使用 scale_y_discrete(labels = paste("Day",mean1 = seq(1,21,length.out =5) ,seq=''))。
ggplot(new_dat, aes(x = Value, y = Class, fill = Class)) + geom_density_ridges() + scale_fill_viridis(discrete = T) + scale_y_discrete(labels = paste("Day",mean1 = seq(1,21,length.out =5) ,seq='')) + theme_bw() + theme(legend.position = "none", panel.grid = element_blank())
此时得到的图形如下：
添加均值线/点
- 第一种方式
使用 ggridges 包中自带的 geom_density_ridges()实现。
ggplot(new_dat, aes(x = Value, y = Class, fill = Class)) + geom_density_ridges() + # 添加每个class的均值 geom_density_ridges(quantile_lines=TRUE, quantile_fun=function(x,...)mean(x))+ scale_fill_viridis(discrete = T) + scale_y_discrete(labels = paste("Day",mean1 = seq(1,21,length.out =5) ,seq='')) + theme_ridges() + theme(legend.position = "none")
- 第二种方式
使用 ggplot2 包中的 stat_summary() 进行添加散点。这里简单形式就以红色点标注，读者也可以根据需求修改颜色的形状大小。
ggplot(new_dat, aes(x = Value, y = Class, fill = Class)) + geom_density_ridges() + stat_summary(fun = "mean", geom = "point", shape = 21, size = 2,fill = "red",color = 'red') + scale_fill_viridis(discrete = T) + scale_y_discrete(labels = paste("Day",mean1 = seq(1,21,length.out =5) ,seq='')) + theme_bw() + theme(legend.position = "none", panel.grid = element_blank())
标注上下分位数
在统计中，上下分位数也是很重要的指标。可以使用 stat_density_ridges() 计算指定分位数 （quantiles = c(0.025, 0.975))）的值并绘制。在通过 scale_fill_manual() 修改填充颜色。
注意：fill 的设置与前面不同，注意区别。
ggplot(new_dat, aes(x = Value, y = Class, fill = factor(after_stat(quantile)))) + stat_density_ridges( geom = "density_ridges_gradient", calc_ecdf = TRUE, quantiles = c(0.025, 0.975)) + scale_fill_manual( name = "Probability", values = c("#FF0000A0", "#A0A0A0A0", "#0000FFA0"), labels = c("(0, 0.025]", "(0.025, 0.975]", "(0.975, 1]") ) + theme_bw() + theme(panel.grid = element_blank())
小编有话说
- 此推文是在小编科研过程中，使用到的图形，可供大家参考。也欢迎大家评论和分享自己的绘图/科研经验。
参考资料
[1]
ggridges: https://cran.r-project.org/web/packages/ggridges/vignettes/introduction.html
