添加p-value和显著性标记：ggsignif和ggpubr
2021-07-09 本文已影响0人
Hayley笔记
在数据分析过程中，常常需要把组间的显著性添加到图形中，但是在ggplot2中实现起来略显麻烦，幸运的是，有很多R包可以帮助我们实现这一操作，比如ggsignif和ggpubr。
1. 统计方法
一般根据数据是否符合正态分布，选择合适的统计方法：
| 统计方法 | 适用情况 |
|---|---|
| t.test() | 比较两组(参数) |
| wilcox.test() | 比较两组(非参数) |
| aov()或anova() | 比较多组(参数) |
| kruskal.test() | 比较多组(非参数) |
2. ggsignif
2.1 ggsignif介绍
ggsignif包主要函数为：geom_signif()和stat_signif()，常用geom_signif()。
# geom_signif参数
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, ...)
| 常用参数 | 说明示范 |
|---|---|
| comparisons | list，设置需要比较的组，比如list(c("a","b"),c("a","c")) |
| test | 上述统计方法，比如t.test() |
| test.args | test传入的参数 |
| map_signif_level | 布尔值，p值直接被当作注释或者以星号替代，比如c(""=0.001,""=0.01,""=0.05) |
| annotations | 带有可选注释的字符向量，如果没有则被忽略 |
| step_increase | 不同组差异标注的距离 |
2.2 ggsignif使用：
安装
# 安装稳定版本：
install.packages("ggsignif")
#安装最新开发版本：
devtools::install_github("const-ae/ggsignif")
使用iris数据集做演示
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的三组两两分别作差异性检验，提前设定好配对分析的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, #修改参数map_signif_level=TRUE
test = wilcox.test)
p1|p2
3. ggpubr
3.1 介绍
ggpubr添加p-value主要使用ggpubr包中的两个函数：compare_means()和stat_compare_mean()。
-
compare_means():可以进行一组或多组间的比较。
compare_means(formula, data, method = "wilcox.test", paired = FALSE,
group.by = NULL, ref.group = NULL, ...)
| 参数 | 含义 |
|---|---|
| formula | 形如x~group，其中x是数值型变量，group是因子，可以是一个或者多个 |
| data | 数据集 |
| method | 比较的方法，默认为"wilcox.test", 其他可选方法为：“t.test”、“anova”、“kruskal.test” |
| paired | 是否要进行paired test(TRUE or FALSE) |
| group_by | 比较时是否要进行分组 |
| ref.group | 是否需要指定参考组 |
-
stat_compare_mean():自动添加p-value、显著性标记到ggplot图中
stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
label = NULL, label.x = NULL, label.y = NULL, ...)
| 参数 | 含义 |
|---|---|
| mapping | 由aes()创建的一套美学映射 |
| comparisons | 指定需要进行比较以及添加p-value、显著性标记的组 |
| hide.ns | 是否要显示显著性标记ns |
| label | 显著性标记的类型，可选项为：p.signif(显著性标记)、p.format(显示p-value) |
| label.x、label.y | 显著性标签调整 |
| … | 其他参数 |
3.2 比较独立的两组
#先加载包
library(ggpubr)
library(patchwork)
data("ToothGrowth") #加载演示数据集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:调整后的p-value。默认为p.adjust.method=“holm”
p.format:四舍五入后的p-value
p.signif:显著性水平
method:用于统计检验的方法
绘制箱线图
p1 <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp",
palette = "lancet", add = "jitter")#添加p-valuep+stat_compare_means()
#使用其他统计检验方法
p2 <- p1+stat_compare_means(method = "t.test")
p1|p2
上述显著性标记可以通过label.x、label.y、hjust及vjust来调整
显著性标记可以通过aes()映射来更改：
aes(label=…p.format…)或aes(lebel=paste0(“p=”,…p.format…))：只显示p-value，不显示统计检验方法
aes(label=…p.signif…)：仅显示显著性水平
aes(label=paste0(…method…,"\n", “p=”,…p.format…))：p-value与显著性水平分行显示
也可以将标签指定为字符向量，不要映射，只需将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,默认方法为”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
#可以通过修改参数label.y来更改标签的位置
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
