数据可视化——R语言为ggplot图形添加P值和显著性水平
本文对一篇英文博客进行翻译，博客原文链接：Add P-values and Significance Levels to ggplots
概述：本文介绍如何轻松地为ggplot图形添加P值和显著性水平：
- 比较两组或多组的均值
- 自动地将P值和显著性水平添加到ggplot图形中，如箱形图，点图，条形图和折线图等
使用工具： R语言中的ggplot2包和ggpubr包
准备
安装和加载R包
本文使用ggpubr包，要求版本高于0.1.3。ggpubr是一个基于ggplot2的计算工具包。
- 直接输入以下命令从CARN中下载安装
install.packages("ggpubr")
- 也可以从GitHub中下载安装最新版本
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
- 加载ggpubr包
library(ggpubr)
ggpubr的官方文档可在以下位置获得：http://www.sthda.com/english/rpkgs/ggpubr
示例数据
示例数据集：ToothGrowth
data("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
均值比较的方法
R中用于两组或多组间均值比较的标准统计方法在之前的文章也有描述：comparing means in R
均值比较的常见方法：
| 方法 | R实现函数 | 描述 |
|---|---|---|
| T-test | t.test() | 比较两组(参数检验) |
| Wilcoxon test | wilcox.test() | 比较两组(非参数检验) |
| ANOVA | aov()或anova() | 比较多组(参数检验) |
| Kruskal-Wallis | kruskal.test() | 比较多组(非参数检验) |
以下链接提供了各种方法的详细介绍：
-
单样本均值与已知均值的比较
-
独立双样本均值比较
-
配对双样本均值比较
-
两组以上的均值比较
- 方差分析（ANOVA, 参数检验）
- Kruskal-Wallis检验（替代单因素方差分析的非参数检验）
用于添加P值的R函数
介绍两个ggpubr包中的函数
- compare_means()：用于执行均值比较
- stat_compare_means()：用于在ggplot图形中自动添加P值和显著性水平
compare_means()
该函数用于执行均值比较。该函数与标准的R函数相比，灵活性更强。
简化形式如下：
compare_means(formula, data, method = "wilcox.test", paired = FALSE,
group.by = NULL, ref.group = NULL, ...)
- formula：指定一个公式，公式形式为 x ~ group，其中，x 表示一个数值型变量，group 表示一个因子型变量，包含一个或多个水平。例如，一个示例公式为 formula = TP53 ~ cancer_group，表示在 cancer_group 对应的各水平间比较TP53的表达水平；也可以同时指定多个响应变量，如 formula = c(TP53, PTEN) ~ cancer_group。
- data：指定一个数据框(data.frame)，数据框需包含formula中的变量。
- method：指定统计检验的方法。默认为“wilcox.test”，即Wilcoxon检验（非参数检验）；也可指定其他统计方法：
- “t.test”，即T检验（参数检验）。“t.test”和“wilcox.test”用于两组样本间的比较。当超过两组时，将会执行两两比较(pairwise comparison)。
- “anova”（参数检验）或 “kruskal.test”（非参数检验），用于执行多组间的单因素方差分析。
- paired：指定一个逻辑变量，表示是否需要执行配对检验，仅适用于t.test 和wilcox.test。
- group.by：指定一个分组变量的字符名，用于在统计检验之前对数据进行分组。当存在group.by指定的变量时，均值比较将在不同水平的各个子集数据中执行。
- ref.group：指定一个组别的字符名，作为对照组(reference group)。如果指定，各个分组水平将与对照组水平进行比较。也可指定ref.group为“.all.”，表示每个分组水平将于所有分组水平（如base-mean）进行比较。
stat_compare_means()
该函数是对ggplot2的扩展，可将均值比较后的P值添加到ggplot图形中，如箱形图、点图、条形图和折线图等。
简化形式如下：
stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
label = NULL, label.x = NULL, label.y = NULL, ...)
-
mapping：通过 aes() 设置绘图时的aesthetic
-
comparisons：指定一个列表(list)，每个列表元素需为长度等于2的向量。向量的内容可以为X轴的两个组别名（字符型），也可以是两个感兴趣组的组别索引（整数值），表示采用指定的两个组别进行比较。
-
hide.ns：逻辑变量，如果设为TRUE，显示显著性水平时将隐藏 ns 字样，即组间差异不显著时不显示 ns 字样。
-
label：指定一个字符串，表示标签类型。可为：“p.signif”（显示显著性水平），“p.format”（显示格式化的P值）。
-
label.x, label.y：指定一个数值，表示显示标签的绝对坐标位置。
-
…：传递给函数compare_means()的参数，如method、paired、ref.group。
独立双样本组间比较
执行统计检验
compare_means(len ~ supp, data = ToothGrowth)
示例结果如下：
## # A tibble: 1 x 8
## .y. group1 group2 p p.adj p.format p.signif method
##
## 1 len OJ VC 0.0645 0.0645 0.064 ns Wilcoxon
method默认为“wilcox.test”（非参数检验），可指定method = “t.test”，表示T检验（参数检验）
返回值为具有以下列的数据框：
-
.y.：用于统计检验的数值变量
-
p：P值
-
p.adj：调整后的P值，调整P值的默认方法为p.adjust.method = “holm”
-
p.format: 格式化的P值
-
p.signif：显著性水平，即用不同数量的 * 表示显著性水平
-
method：用于组间比较的统计方法
创建添加P值的箱形图
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
color = "supp", palette = "jco",
add = "jitter")
# Add p-value
p + stat_compare_means()
# Change method
p + stat_compare_means(method = "t.test")
注意：显示P值的标签位置可以通过如下参数来调整：label.x, label.y, hjust 和vjust。
显示P值的标签默认为 compare_means() 返回值中的 method 和 p 的组合。也可以通过 aes() 函数指定为其他显示形式。例如：
aes(label = ..p.format..) 或 aes(label = paste0(“p =”, ..p.format..))
表示只显示格式化的P值，而不显示method
aes(label = ..p.signif..)
表示仅显示显著性水平
aes(label = paste0(..method.., “\n”, “p =”, ..p.format..))
表示在method名和P值之间添加换行符(“\n”)
示例如下：
p + stat_compare_means( aes(label = ..p.signif..),
label.x = 1.5, label.y = 40)
另外，也可以将参数label指定为字符向量：
p + stat_compare_means( label = "p.signif", label.x = 1.5, label.y = 40)
配对双样本组间比较
执行统计检验
compare_means(len ~ supp, data = ToothGrowth, paired = TRUE)
示例结果如下：
## # A tibble: 1 x 8
## .y. group1 group2 p p.adj p.format p.signif method
##
## 1 len OJ VC 0.00431 0.00431 0.0043 ** Wilcoxon
使用函数 ggpaired() 可视化配对数据
ggpaired(ToothGrowth, x = "supp", y = "len",
color = "supp", line.color = "gray", line.size = 0.4,
palette = "jco")+
stat_compare_means(paired = TRUE)
多组样本的组间比较
- 全局检验（所有组的均值比较）
# Global test
compare_means(len ~ dose, data = ToothGrowth, method = "anova")
示例结果如下：
## # A tibble: 1 x 6
## .y. p p.adj p.format p.signif method
##
## 1 len 9.53e-16 9.53e-16 9.5e-16 **** Anova
添加全局检验的P值（所有组比较总的P值）
# Default method = "kruskal.test" for multiple groups
ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means()
# Change method to anova
ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(method = "anova")
- 两两比较(Pairwise comparisons)
如果分组变量包含两个以上的水平，两两比较的检验(pairwise test)将自动执行。默认方法为“wilcox.test”，也可设置为“t.test”。
# Perorm pairwise comparisons
compare_means(len ~ dose, data = ToothGrowth)
示例结果如下：
## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
##
## 1 len 0.5 1 7.02e-06 1.40e-05 7.0e-06 **** Wilcoxon
## 2 len 0.5 2 8.41e-08 2.52e-07 8.4e-08 **** Wilcoxon
## 3 len 1 2 1.77e-04 1.77e-04 0.00018 *** Wilcoxon
# Visualize: Specify the comparisons you want
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
stat_compare_means(label.y = 50) # Add global p-value
如果需要指定标签显示的Y轴位置，可使用参数label.y
ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(comparisons = my_comparisons, label.y = c(29, 35, 40))+
stat_compare_means(label.y = 45)
注：ggsignif 包也可以很方便的为条形图添加组间比较的P值
- 相对于对照组的多重两两比较的检验(Multiple pairwise tests)
# Pairwise comparison against reference
compare_means(len ~ dose, data = ToothGrowth, ref.group = "0.5",
method = "t.test")
示例结果如下：
## # A tibble: 2 x 8
## .y. group1 group2 p p.adj p.format p.signif method
##
## 1 len 0.5 1 6.70e-09 6.70e-09 6.7e-09 **** T-test
## 2 len 0.5 2 1.47e-16 2.94e-16 < 2e-16 **** T-test
# Visualize
ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(method = "anova", label.y = 40)+ # Add global p-value
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = "0.5") # Pairwise comparison against reference
- 相对于所有组(base-mean)的多重两两比较的检验
# Comparison of each group against base-mean
compare_means(len ~ dose, data = ToothGrowth, ref.group = ".all.",
method = "t.test")
示例结果如下：
## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
##
## 1 len .all. 0.5 1.24e-06 3.73e-06 1.2e-06 **** T-test
## 2 len .all. 1 5.67e-01 5.67e-01 0.57 ns T-test
## 3 len .all. 2 1.37e-05 2.74e-05 1.4e-05 **** T-test
# Visualize
ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(method = "anova", label.y = 40)+ # Add global p-value
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = ".all.") # Pairwise comparison against all
下面使用Github中可用的骨髓瘤数据集展示一些典型情况，其中，与“.all.”的比较将很有用。
将依据患者molecular进行分组，绘制各个组别DEPDC1基因的表达水平。目的是比较各个组别之间是否存在差别，如果有差别，差别又在哪里。
要回答以上问题，可以在所有7个组之间进行两两比较(pairwise comparison)。 由于组别较多，将导致很多种组别组合，这将很难解释。
一个简单的解决办法是将7组中的每一组与“.all.”（如base-mean）比较。当检验结果显著时，可以得出结论：与所有组相比，xxx组的DEPDC1的表达水平显著降低或显著升高。
# Load myeloma data from GitHub
myeloma <- read.delim("https://raw.githubusercontent.com/kassambara/data/master/myeloma.txt")
# Perform the test
compare_means(DEPDC1 ~ molecular_group, data = myeloma,
ref.group = ".all.", method = "t.test")
示例结果如下：
## # A tibble: 7 x 8
## .y. group1 group2 p p.adj p.format p.signif method
##
## 1 DEPDC1 .all. Cyclin D-1 0.149690 0.44907 0.14969 ns T-test
## 2 DEPDC1 .all. Cyclin D-2 0.523143 1.00000 0.52314 ns T-test
## 3 DEPDC1 .all. Hyperdiploid 0.000282 0.00169 0.00028 *** T-test
## 4 DEPDC1 .all. Low bone disease 0.005084 0.02542 0.00508 ** T-test
## 5 DEPDC1 .all. MAF 0.086107 0.34443 0.08611 ns T-test
## 6 DEPDC1 .all. MMSET 0.576291 1.00000 0.57629 ns T-test
## # ... with 1 more rows
注意：上面的R代码中myeloma数据的下载地址存在错误，可从以下链接找到myeloma数据：https://github.com/kassambara/data/blob/master/myeloma.txt
然后自己将全部数据复制并保存为myeloma.txt，从本地读取myeloma.txt。
# Visualize the expression profile
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1", color = "molecular_group",
add = "jitter", legend = "none") +
rotate_x_text(angle = 45)+
geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2)+ # Add horizontal line at base mean
stat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-value
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = ".all.") # Pairwise comparison against all
根据上图的结果，可以得出proliferation组的DEPDC1表达水平显著升高；Hyperdiploid组和Low bone disease组中DEPDC1表达水平显著降低。
注意：想要隐藏 ns 标志，可以设置参数 hide.ns = TRUE
# Visualize the expression profile
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1", color = "molecular_group",
add = "jitter", legend = "none") +
rotate_x_text(angle = 45)+
geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2)+ # Add horizontal line at base mean
stat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-value
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = ".all.", hide.ns = TRUE) # Pairwise comparison against all
多个分组变量
- 使用另一个变量进行分组后再执行独立双样本比较
执行统计检验：
compare_means(len ~ supp, data = ToothGrowth,
group.by = "dose")
示例结果如下：
## # A tibble: 3 x 9
## dose .y. group1 group2 p p.adj p.format p.signif method
##
## 1 0.5 len OJ VC 0.02319 0.0464 0.023 * Wilcoxon
## 2 1.0 len OJ VC 0.00403 0.0121 0.004 ** Wilcoxon
## 3 2.0 len OJ VC 1.00000 1.0000 1.000 ns Wilcoxon
在上面的示例中，对于分类变量dose的每一个水平，分类变量supp又将数据分为两个水平：OJ和VC，然后在这两个水平上对数值变量len进行均值比较。
可视化：创建一个按组划分的多面板框图(此处为“dose”)
# Box plot facetted by "dose"
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
color = "supp", palette = "jco",
add = "jitter",
facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format")
# Or use significance symbol as label
p + stat_compare_means(label = "p.signif", label.x = 1.5)
注意：想要隐藏 ns 标志，可以设置参数hide.ns = TRUE
可视化：创建一个包含所有箱形图的单一面板。X表示dose，Y表示len，颜色表示supp
p <- ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "supp", palette = "jco",
add = "jitter")
p + stat_compare_means(aes(group = supp))
# Show only p-value
p + stat_compare_means(aes(group = supp), label = "p.format")
# Use significance symbol as label
p + stat_compare_means(aes(group = supp), label = "p.signif")
- 使用另一个变量进行分组后再执行配对双样本比较
执行统计检验：
compare_means(len ~ supp, data = ToothGrowth,
group.by = "dose", paired = TRUE)
实例结果如下：
## # A tibble: 3 x 9
## dose .y. group1 group2 p p.adj p.format p.signif method
##
## 1 0.5 len OJ VC 0.0330 0.0659 0.033 * Wilcoxon
## 2 1.0 len OJ VC 0.0191 0.0572 0.019 * Wilcoxon
## 3 2.0 len OJ VC 1.0000 1.0000 1.000 ns Wilcoxon
可视化：创建一个按组划分的多面板框图(此处为“dose”)
# Box plot facetted by "dose"
p <- ggpaired(ToothGrowth, x = "supp", y = "len",
color = "supp", palette = "jco",
line.color = "gray", line.size = 0.4,
facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format", paired = TRUE)
其他绘图方式
- 条形图和折线图（一个分组变量）
# Bar plot of mean +/-se
ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
stat_compare_means() + # Global p-value
stat_compare_means(ref.group = "0.5", label = "p.signif",
label.y = c(22, 29)) # compare to ref.group
# Line plot of mean +/-se
ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
stat_compare_means() + # Global p-value
stat_compare_means(ref.group = "0.5", label = "p.signif",
label.y = c(22, 29))
- 条形图和折线图（两个分组变量）
ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se",
color = "supp", palette = "jco",
position = position_dodge(0.8))+
stat_compare_means(aes(group = supp), label = "p.signif", label.y = 29)
ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se",
color = "supp", palette = "jco")+
stat_compare_means(aes(group = supp), label = "p.signif",
label.y = c(16, 25, 29))
注意：经过实际测试，笔者发现R语言中的统计方法计算结果的P值与SPSS中的P值存在差异。如，常规的方差分析(ANOVA) + 事后两两组间比较（如Bonferroni校正）使用上述R函数就很难得出与SPSS中一致的结果。如果需要使用SPSS的统计P值，建议对生成的图形进行后期修改。
References
扫一扫
4707
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
