基于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 天的后验结果。这里的数据不算是“整洁 Tidy”数据。需要通过 pivot_longer() 进行数据整理。

dat %>% pivot_longer(`z1`:`z5`,
                     names_to = "Class",
                     values_to = "Value"
                     ) -> new_dat

我们将 z1:z5 名称存储到 "Class" 中，数值存储到 "Value" 中。转化后的数据如下所示：

初级版本

根据以上数据 new_dat，我们先使用 ggridges 包中的 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())

小编有话说

• 此推文是在小编科研过程中，使用到的图形，可供大家参考。也欢迎大家评论和分享自己的绘图/科研经验。