在數(shù)據(jù)分析過程中，常常需要把組間的顯著性添加到圖形中，但是在ggplot2中實現(xiàn)起來略顯麻煩，幸運的是，有很多R包可以幫助我們實現(xiàn)這一操作，比如ggsignif和ggpubr。
1. 統(tǒng)計方法
一般根據(jù)數(shù)據(jù)是否符合正態(tài)分布，選擇合適的統(tǒng)計方法：
| 統(tǒng)計方法 | 適用情況 |
|---|---|
| t.test() | 比較兩組(參數(shù)) |
| wilcox.test() | 比較兩組(非參數(shù)) |
| aov()或anova() | 比較多組(參數(shù)) |
| kruskal.test() | 比較多組(非參數(shù)) |
2. ggsignif
2.1 ggsignif介紹
ggsignif包主要函數(shù)為：geom_signif()和stat_signif()，常用geom_signif()。
# geom_signif參數(shù)
geom_signif(mapping = NULL, data = NULL, stat = "signif",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, comparisons = NULL, test = "wilcox.test",
test.args = NULL, annotations = NULL, map_signif_level = FALSE,
y_position = NULL, xmin = NULL, xmax = NULL, margin_top = 0.05,
step_increase = 0, tip_length = 0.03, size = 0.5, textsize = 3.88,
family = "", vjust = 0, ...)
| 常用參數(shù) | 說明示范 |
|---|---|
| comparisons | list，設(shè)置需要比較的組，比如list(c("a","b"),c("a","c")) |
| test | 上述統(tǒng)計方法，比如t.test() |
| test.args | test傳入的參數(shù) |
| map_signif_level | 布爾值，p值直接被當(dāng)作注釋或者以星號替代，比如c(""=0.001,""=0.01,""=0.05) |
| annotations | 帶有可選注釋的字符向量，如果沒有則被忽略 |
| step_increase | 不同組差異標(biāo)注的距離 |
2.2 ggsignif使用：
安裝
# 安裝穩(wěn)定版本：
install.packages("ggsignif")
#安裝最新開發(fā)版本：
devtools::install_github("const-ae/ggsignif")
使用iris數(shù)據(jù)集做演示
library(ggplot2)
library(ggthemr) #載入主題配置包
library(ggsignif) #載入ggsignif
library(patchwork) #載入拼圖包
head(iris)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species
# 1 5.1 3.5 1.4 0.2 setosa
# 2 4.9 3.0 1.4 0.2 setosa
# 3 4.7 3.2 1.3 0.2 setosa
# 4 4.6 3.1 1.5 0.2 setosa
# 5 5.0 3.6 1.4 0.2 setosa
# 6 5.4 3.9 1.7 0.4 setosa
Species的三組兩兩分別作差異性檢驗，提前設(shè)定好配對分析的list
compaired <- list(c("versicolor", "virginica"),
c("versicolor","setosa"),
c("virginica","setosa"))
繪制geom_boxplot()和小提琴圖geom_violin()
ggthemr("flat")
p1 <- ggplot(iris, aes(Species, Sepal.Width, fill = Species)) +
geom_boxplot() +
ylim(1.5, 6.5) +
geom_signif(comparisons = compaired,
step_increase = 0.3,
map_signif_level = F,
test = wilcox.test)
p2 <- ggplot(iris, aes(Species, Sepal.Width, fill = Species)) +
geom_violin() +
ylim(1.5, 6.5) +
geom_signif(comparisons = compaired,
step_increase = 0.3,
map_signif_level = T, #修改參數(shù)map_signif_level=TRUE
test = wilcox.test)
p1|p2
3. ggpubr
3.1 介紹
ggpubr添加p-value主要使用ggpubr包中的兩個函數(shù)：compare_means()和stat_compare_mean()。
-
compare_means():可以進行一組或多組間的比較。
compare_means(formula, data, method = "wilcox.test", paired = FALSE,
group.by = NULL, ref.group = NULL, ...)
| 參數(shù) | 含義 |
|---|---|
| formula | 形如x~group，其中x是數(shù)值型變量，group是因子，可以是一個或者多個 |
| data | 數(shù)據(jù)集 |
| method | 比較的方法，默認(rèn)為"wilcox.test", 其他可選方法為：“t.test”、“anova”、“kruskal.test” |
| paired | 是否要進行paired test(TRUE or FALSE) |
| group_by | 比較時是否要進行分組 |
| ref.group | 是否需要指定參考組 |
-
stat_compare_mean():自動添加p-value、顯著性標(biāo)記到ggplot圖中
stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
label = NULL, label.x = NULL, label.y = NULL, ...)
| 參數(shù) | 含義 |
|---|---|
| mapping | 由aes()創(chuàng)建的一套美學(xué)映射 |
| comparisons | 指定需要進行比較以及添加p-value、顯著性標(biāo)記的組 |
| hide.ns | 是否要顯示顯著性標(biāo)記ns |
| label | 顯著性標(biāo)記的類型，可選項為：p.signif(顯著性標(biāo)記)、p.format(顯示p-value) |
| label.x、label.y | 顯著性標(biāo)簽調(diào)整 |
| … | 其他參數(shù) |
3.2 比較獨立的兩組
#先加載包
library(ggpubr)
library(patchwork)
data("ToothGrowth") #加載演示數(shù)據(jù)集ToothGrowth
head(ToothGrowth)
# len supp dose
# 1 4.2 VC 0.5
# 2 11.5 VC 0.5
# 3 7.3 VC 0.5
# 4 5.8 VC 0.5
# 5 6.4 VC 0.5
# 6 10.0 VC 0.5
compare_means(len~supp, data=ToothGrowth)
## A tibble: 1 x 8
# .y. group1 group2 p p.adj p.format p.signif method
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
#1 len OJ VC 0.0645 0.064 0.064 ns Wilcox…
y:測試中使用的y變量
p:p-value
p.adj:調(diào)整后的p-value。默認(rèn)為p.adjust.method=“holm”
p.format:四舍五入后的p-value
p.signif:顯著性水平
method:用于統(tǒng)計檢驗的方法
繪制箱線圖
p1 <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp",
palette = "lancet", add = "jitter")#添加p-valuep+stat_compare_means()
#使用其他統(tǒng)計檢驗方法
p2 <- p1+stat_compare_means(method = "t.test")
p1|p2
上述顯著性標(biāo)記可以通過label.x、label.y、hjust及vjust來調(diào)整
顯著性標(biāo)記可以通過aes()映射來更改：
aes(label=…p.format…)或aes(lebel=paste0(“p=”,…p.format…))：只顯示p-value，不顯示統(tǒng)計檢驗方法
aes(label=…p.signif…)：僅顯示顯著性水平
aes(label=paste0(…method…,"\n", “p=”,…p.format…))：p-value與顯著性水平分行顯示
也可以將標(biāo)簽指定為字符向量，不要映射，只需將p.signif兩端的…去掉即可
3.3 配對樣本
compare_means(len~supp, data=ToothGrowth, paired = TRUE)
## A tibble: 1 x 8
# .y. group1 group2 p p.adj p.format p.signif
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
#1 len OJ VC 0.00431 0.0043 0.0043 **
## … with 1 more variable: method <chr>
利用ggpaired()進行可視化
ggpaired(ToothGrowth, x="supp", y="len", color = "supp", line.color = "gray",
line.size = 0.4, palette = "jco")+ stat_compare_means(paired = TRUE)
3.4 多組比較 Global test
compare_means(len~dose, data=ToothGrowth, method = "anova")
## A tibble: 1 x 6
# .y. p p.adj p.format p.signif method
# <chr> <dbl> <dbl> <chr> <chr> <chr>
#1 len 9.53e-16 9.5e-16 9.5e-16 **** Anova
繪圖
p1=ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+
stat_compare_means()
#使用其他的方法
p2=ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+
stat_compare_means(method = "anova")
p1|p2
Pairwise comparisons:如果分組變量中包含兩個以上的水平，那么會自動進行pairwise test,默認(rèn)方法為”wilcox.test”
compare_means(len~dose, data=ToothGrowth)
## A tibble: 3 x 8
# .y. group1 group2 p p.adj p.format p.signif
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
#1 len 0.5 1 7.02e-6 1.4e-5 7.0e-06 ****
#2 len 0.5 2 8.41e-8 2.5e-7 8.4e-08 ****
#3 len 1 2 1.77e-4 1.8e-4 0.00018 ***
## … with 1 more variable: method <chr>
#可以指定比較哪些組
my_comparisons <- list(c("0.5", "1"), c("1", "2"), c("0.5", "2"))
p1=ggboxplot(ToothGrowth, x="dose", y="len", color = "dose",palette = "jco")+
stat_compare_means(comparisons=my_comparisons)+ # Add pairwisecomparisons p-value
stat_compare_means(label.y = 50) # Add global p-value
#可以通過修改參數(shù)label.y來更改標(biāo)簽的位置
p2=ggboxplot(ToothGrowth, x="dose", y="len", color = "dose",palette = "jco")+
stat_compare_means(comparisons=my_comparisons, label.y = c(29, 35, 40))+
# Add pairwise comparisons p-value
stat_compare_means(label.y = 45) # Add global p-value
p1|p2
