優(yōu)雅地復現(xiàn)R統(tǒng)計圖：從「梵婀玲」圖（violin plot）開始
作者：狐貝貝
? ? ? ? ANOVA（方差分析）或Kruskal-Wallis H秩和檢驗是統(tǒng)計學中用于檢驗兩個或多個樣本均值/中位數(shù)之間是否存在顯著差異的常用算法。
? ? ? ? 在統(tǒng)計學的學習資料或論文當中，我們經(jīng)常能看到如下圖所示ANOVA分析或Kruskal-Wallis秩和檢驗的小提琴圖。
? ? ? ? ?小提琴圖清楚地展示了各組的數(shù)據(jù)分布，ANOVA分析或Kruskal-Wallis秩和檢驗的結果，以及兩兩組間比較的顯著性。
對于復現(xiàn)這樣一張簡單的統(tǒng)計圖，同學們是不是摩拳擦掌，躍躍欲試呢？
? ? ? ?下面，小猿團隊的小狐帶您一步步使用R語言的ggplot函數(shù)或統(tǒng)計猿在線工具實現(xiàn)Kruskal-Wallis秩和檢驗及數(shù)據(jù)的可視化。本講內(nèi)容包括數(shù)據(jù)準備與整理、Kruskal-Wallis秩和檢驗和事后多重比較檢驗。
?在學習教程之前，請同學們自行完成R語言GUI和RStudio兩個軟件的安裝。
未安裝“ggplot2”等R包的同學請先通過install.packages("包的名字")自行安裝。
Step 1:構造模擬數(shù)據(jù)
首先我們用runif函數(shù)構造了5組模擬數(shù)據(jù)，每組10個隨機數(shù)。runif函數(shù)可以生成指定范圍內(nèi)的均勻分布隨機數(shù)。
library(tidyverse)
library(ggpubr)
library(ggsignif)
library(rstatix)
#隨機數(shù)生成每次會有所不同
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:寬數(shù)據(jù)轉變?yōu)殚L數(shù)據(jù)，以及保存數(shù)據(jù)
通過pivot_longer函數(shù)將寬格式的數(shù)據(jù)轉換為長格式。長格式更方便進行分組分析。
由于代碼生成的數(shù)據(jù)具有隨機性，每次都不一樣，同學們可以保存數(shù)據(jù)，下次直接讀取數(shù)據(jù)。
data_long <- pivot_longer(data, cols = everything(),
????????????????????????? names_to = "group", values_to = "Score") %>%
? as.data.frame()
#保存數(shù)據(jù)
write.csv(data_long,"data_long.csv")
?
#讀取數(shù)據(jù)
data_long=read.csv("data_long.csv")
#數(shù)據(jù)進行因子化
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", #根據(jù)校正后的p value換算顯著性符號
#對于顯著性符號，小編一般只用到以下幾個層級（做多3個星號，默認最多是4個星號）
?????????????????? cutpoints = c(0, 0.001, 0.01, 0.05, 1),
?????????????????? symbols = c("***", "**", "*", "ns"))
??
Step 5：繪制小提琴圖及添加顯著性符號
核密度圖、散點圖、箱線圖共同構成了小提琴圖，直觀展示了不同組的分布情況。
此外，添加了統(tǒng)計檢驗的顯著性標記。
?
#使用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",
??? #去掉網(wǎng)格線
??? 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
# 根據(jù)顯著性檢驗結果，添加顯著性標記：
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循環(huán)添加了顯著性符號，完成小提琴圖的繪制。
使用for循環(huán)算法，完成顯著性符號的添加。通過顯著性符號的顏色，即可判斷出是與哪一組數(shù)據(jù)進行比較。
例如，TIME_IA列上方的顯著性符號*為靛青色，說明經(jīng)校正假陽性的dunn’s檢驗反映出，TIME_IA組顯著高于（高低是主要是基于中位數(shù)的判斷）TIME_IE組。
?
Step 6: 將圖片保存文件
保存為矢量圖格式，后續(xù)可進行字體、字號、線條粗細等細節(jié)的調(diào)整。
ggsave("single_plot.pdf", height = 4, width = 4)
ggsave("single_plot.svg", height = 4, width = 4)
通過顯著性符號的顏色，即可判斷出是與哪一組數(shù)據(jù)進行比較。
例如，TIME_IA列上方的顯著性符號*為靛青色，說明經(jīng)校正假陽性的dunn’s檢驗反映出，TIME_IA組顯著高于（高低是主要是基于中位數(shù)的判斷）TIME_IE組。
?
Step 7: 將圖片保存文件
保存為矢量圖格式，后續(xù)可進行字體、字號、線條粗心等細節(jié)的調(diào)整。
ggsave("single_plot.pdf", height = 4, width = 4)
ggsave("single_plot.svg", height = 4, width = 4)
如果認為寫代碼比較困難，我們還可以再嘗試用統(tǒng)計猿軟件實現(xiàn)統(tǒng)計圖的復現(xiàn)。
打開我們的統(tǒng)計猿網(wǎng)站（https://fast.statsape.com/）。
選擇高級制圖。
點開小提琴圖-多重比較小提琴圖。
選擇上傳本地數(shù)據(jù)。還未生成圖片之前，左側看到的小提琴圖是樣例展示，非真實數(shù)據(jù)（可觀察橫坐標）。
點擊文件夾圖標進行上傳數(shù)據(jù)表（剛剛使用R語言生成的數(shù)據(jù)）。提示“上傳成功”。
在右側設置好各項屬性后，包括橫坐標、縱坐標、分類變量、配色等。小編作圖時，各項參數(shù)均與R語言中的參數(shù)保持了一致。
點擊“開始”，不一會兒便生成了圖片。
在本篇教程中，我們實現(xiàn)了Kruskal-Wallis秩和檢驗分析的完整流程,包括假設檢驗、可視化展示,對學習統(tǒng)計分析方法很有幫助。同學們在實踐中，可替換為自己的數(shù)據(jù)，進行統(tǒng)計學分析。
