基于 ggridges 绘制剩余使用寿命密度图
点击下方公众号，回复资料分享，收获惊喜
简介
最近科研过程中，需要可视化系统剩余使用寿命（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())
小编有话说
此推文是在小编科研过程中，使用到的图形，可供大家参考。也欢迎大家评论和分享自己的绘图/科研经验。
参考资料
ggridges: https://cran.r-project.org/web/packages/ggridges/vignettes/introduction.html
推荐： 可以保存以下照片，在b站扫该二维码，或者b站搜索【庄闪闪】观看Rmarkdown系列的视频教程。Rmarkdown视频新增两节视频（写轮眼幻灯片制作）需要视频内的文档，可在公众号回复【rmarkdown】
R沟通｜Rmarkdown教程（4）
R沟通｜Rmarkdown教程（3）
R沟通｜Rmarkdown教程（2）