优雅地复现R统计图：从「梵婀玲」图（violin plot）开始

作者：狐贝贝

        ANOVA（方差分析）或Kruskal-Wallis H秩和检验是统计学中用于检验两个或多个样本均值/中位数之间是否存在显著差异的常用算法。

        在统计学的学习资料或论文当中，我们经常能看到如下图所示ANOVA分析或Kruskal-Wallis秩和检验的小提琴图。

         小提琴图清楚地展示了各组的数据分布，ANOVA分析或Kruskal-Wallis秩和检验的结果，以及两两组间比较的显著性。

对于复现这样一张简单的统计图，同学们是不是摩拳擦掌，跃跃欲试呢？

       下面，小猿团队的小狐带您一步步使用R语言的ggplot函数或统计猿在线工具实现Kruskal-Wallis秩和检验及数据的可视化。本讲内容包括数据准备与整理、Kruskal-Wallis秩和检验和事后多重比较检验。

 在学习教程之前，请同学们自行完成R语言GUI和RStudio两个软件的安装。

未安装“ggplot2”等R包的同学请先通过install.packages("包的名字")自行安装。

Step 1:构造模拟数据

首先我们用runif函数构造了5组模拟数据，每组10个随机数。runif函数可以生成指定范围内的均匀分布随机数。

library(tidyverse)

library(ggpubr)

library(ggsignif)

library(rstatix)

#随机数生成每次会有所不同

data <- data.frame(

  TIME_IA = runif(10, min = 0.05, max = 0.4),

  TIME_ISM = runif(10, min = -0.2, max = 0.1),

  TIME_ISS = runif(10, min = 0.1, max = 0.5),

  TIME_IE = runif(10, min = -0.25, max = 0.2),

  TIME_IR = runif(10, min = 0.2, max = 0.6)

)

Step 2:宽数据转变为长数据，以及保存数据

通过pivot_longer函数将宽格式的数据转换为长格式。长格式更方便进行分组分析。

由于代码生成的数据具有随机性，每次都不一样，同学们可以保存数据，下次直接读取数据。

data_long <- pivot_longer(data, cols = everything(),

                          names_to = "group", values_to = "Score") %>%

  as.data.frame()

#保存数据

write.csv(data_long,"data_long.csv")

 

#读取数据

data_long=read.csv("data_long.csv")

#数据进行因子化

data_long$group <- factor(data_long$group, levels = colnames(data))

 

Step 3：Kruskal-Wallis秩和检验

 

使用kruskal_test进行Kruskal-Wallis秩和检验。p值小于0.05，说明至少有两组之间存在显著差异，但具体是哪两组需要进行事后多重比较。

kw_test <- data_long %>%

  kruskal_test(Score ~ group)

kw_test

Step 4：事后多重比较检验

采用dunn_test（Dunn’s test）进行事后多重比较,并用bonferroni方法控制整体型I错误。标注出了不同组间的显著性差异。

stat.test <- data_long %>%

  dunn_test(

    Score ~ group,

p.adjust.method = "Bonferroni

"

  )%>%

  add_significance("p.adj", #根据校正后的p value换算显著性符号

#对于显著性符号，小编一般只用到以下几个层级（做多3个星号，默认最多是4个星号）

                   cutpoints = c(0, 0.001, 0.01, 0.05, 1),

                   symbols = c("***", "**", "*", "ns"))

  

Step 5：绘制小提琴图及添加显著性符号

核密度图、散点图、箱线图共同构成了小提琴图，直观展示了不同组的分布情况。

此外，添加了统计检验的显著性标记。

 

#使用ggsci包生成5个《自然》期刊颜色

library("scales")

library(ggsci)

 

show_col(pal_npg("nrc")(nlevels(data_long$group)))

 

colors=pal_npg("nrc")(nlevels(data_long$group))

colors_a=pal_npg(palette = c("nrc"), alpha = 0.6)(nlevels(data_long$group))

 

#写作图的主题

#喜欢拐角（GraphPad Prism）风格的主题用这个

prism=  theme_classic()+

  theme(

    legend.position = "none",

    #tick深灰色

    axis.ticks = element_line(color="grey20"),

    axis.title.x=element_text(size = 8,color="black",face="bold"),

    axis.title.y=element_text(size = 8,color="black",face="bold"),

    legend.title = element_text(size = 8,color="black",face="bold"),

    axis.text.x=element_text(

      angle = 45, #横坐标旋转45°

      hjust = 1,

      vjust=1,

      size = 8,color="black",

      margin=unit(c(0.2,0.2,0.2,0.2), "cm")),

    axis.text.y = element_text(

     

      size = 8,color="black",

      margin=unit(c(0.2,0.2,0.2,0.2), "cm")),

    legend.text = element_text(size = 8,color="black"),

    axis.ticks.length=unit(0.2, "cm")#tick长度

  )

 

#喜欢完整包围起来风格的主题用这个

bw=theme_bw()+

  theme(

    legend.position = "none",

    #去掉网格线

    panel.grid = element_blank(),

    #tick深灰色

    axis.ticks = element_line(color="grey20"),

    axis.title.x=element_text(size = 8,color="black",face="bold"),

    axis.title.y=element_text(size = 8,color="black",face="bold"),

    legend.title = element_text(size = 8,color="black",face="bold"),

    axis.text.x=element_text(

      angle = 45, #横坐标旋转45°

      hjust = 1,

      vjust=1,

      size = 8,color="black",

      margin=unit(c(0.2,0.2,0.2,0.2), "cm")),

    axis.text.y = element_text(

     

      size = 8,color="black",

      margin=unit(c(0.2,0.2,0.2,0.2), "cm")),

    legend.text = element_text(size = 8,color="black"),

    axis.ticks.length=unit(0.2, "cm")#tick长度

  )

 

#先绘制图像

p <- ggplot()+

  #核密度图：

  geom_violin(data=data_long,aes(group, Score, color = group, fill = group), width = 0.75)+

    # 箱线图：

  geom_boxplot(data=data_long,aes(group, Score, color = group), width = 0.2)+

  # 抖动散点：

  geom_jitter(data=data_long,aes(group, Score, color = group), width = 0.05)+

 

  #添加KW检验结果

  geom_text(aes(x=1.5,y=0.6),label=paste0(kw_test$method,", ",

                                        ifelse(kw_test$p<0.001,"p < 0.001",paste0("p = ", sprintf("%.4f", kw_test$p)))))+

  # 颜色模式：

  scale_color_manual(values = colors )+

  scale_fill_manual(values = colors_a )+

    scale_y_continuous(limits = c(-0.3,0.62),breaks=seq(-0.3,0.6,0.15))+

  xlab("Group")+

  # 主题

  bw

p

# 根据显著性检验结果，添加显著性标记：

x_value <- rep(1:4, 4:1)

y_value <- rep(apply(data, 2, max)[1:4], 4:1) + 0.01

y_value <- y_value + c(0.03*1:4, 0.03*1:3, 0.03*1:2, 0.03)

color_value <- c(colors[2:5], colors[3:5], colors[4:5], colors[5])

 

for (i in 1:nrow(stat.test)) {

  if (stat.test$p.adj.signif[i] != "ns") {

    y_tmp <- y_value[i]

    p <- p+annotate(geom = "text",

                    label = stat.test$p.adj.signif[i],

                    x = x_value[i],

                    y = y_tmp,

                    color = color_value[i])

  }

}

 

p

通过for循环添加了显著性符号，完成小提琴图的绘制。

使用for循环算法，完成显著性符号的添加。通过显著性符号的颜色，即可判断出是与哪一组数据进行比较。

例如，TIME_IA列上方的显著性符号*为靛青色，说明经校正假阳性的dunn’s检验反映出，TIME_IA组显著高于（高低是主要是基于中位数的判断）TIME_IE组。

 

Step 6: 将图片保存文件

保存为矢量图格式，后续可进行字体、字号、线条粗细等细节的调整。

ggsave("single_plot.pdf", height = 4, width = 4)

ggsave("single_plot.svg", height = 4, width = 4)

通过显著性符号的颜色，即可判断出是与哪一组数据进行比较。

例如，TIME_IA列上方的显著性符号*为靛青色，说明经校正假阳性的dunn’s检验反映出，TIME_IA组显著高于（高低是主要是基于中位数的判断）TIME_IE组。

 

Step 7: 将图片保存文件

保存为矢量图格式，后续可进行字体、字号、线条粗心等细节的调整。

ggsave("single_plot.pdf", height = 4, width = 4)

ggsave("single_plot.svg", height = 4, width = 4)

如果认为写代码比较困难，我们还可以再尝试用统计猿软件实现统计图的复现。

打开我们的统计猿网站（https://fast.statsape.com/）。

选择高级制图。

点开小提琴图-多重比较小提琴图。

选择上传本地数据。还未生成图片之前，左侧看到的小提琴图是样例展示，非真实数据（可观察横坐标）。

点击文件夹图标进行上传数据表（刚刚使用R语言生成的数据）。提示“上传成功”。

在右侧设置好各项属性后，包括横坐标、纵坐标、分类变量、配色等。小编作图时，各项参数均与R语言中的参数保持了一致。

点击“开始”，不一会儿便生成了图片。

在本篇教程中，我们实现了Kruskal-Wallis秩和检验分析的完整流程,包括假设检验、可视化展示,对学习统计分析方法很有帮助。同学们在实践中，可替换为自己的数据，进行统计学分析。