R语言添加p-value和显著性标记

2024-02-01 01:40
文章标签 语言 标记 value 显著性

本文主要是介绍R语言添加p-value和显著性标记，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

R语言添加p-value和显著性标记，原文链接 https://mp.weixin.qq.com/s/gRw0krS3LY7c0QK9y47EJw

作者简介 Introduction

taoyan：伪码农，R语言爱好者，爱开源。

个人博客: https://ytlogos.github.io/

往期回顾

  • 一条命令轻松绘制CNS顶级配图-ggpubr
  • R语言聚类分析–cluster, factoextra

上篇文章中提了一下如何通过ggpubr包为ggplot图添加p-value以及显著性标记,本文将详细介绍。利用数据集ToothGrowth进行演示。

ggpubr安装和加载

# 直接从CRAN安装
install.packages("ggpubr", repo="http://cran.us.r-project.org")#先加载包
library(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中常用的比较方法主要有下面几种：

方法R函数描述
T-testt.test()比较两组(参数)
Wilcoxon testwilcox.test()比较两组(非参数)
ANOVAaov()或anova()比较多组(参数)
Kruskal-Walliskruskal.test()比较多组(非参数)

添加p-value

主要利用ggpubr包中的两个函数：

  • compare_means():可以进行一组或多组间的比较

  • stat_compare_mean():自动添加p-value、显著性标记到ggplot图中

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_means()函数

主要用法：

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:显著性标签调整

…:其他参数

比较独立的两组

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.06449067 0.06449067    0.064       ns Wilcoxon

结果解释：

.y:测试中使用的y变量

p:p-value

p.adj:调整后的p-value。默认为p.adjust.method=”holm”

p.format:四舍五入后的p-value

p.signif:显著性水平

method:用于统计检验的方法

绘制箱线图

# 绘制箱线图
p <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp",palette = "jco", add = "jitter") 
# 添加p-value, 默认是Wilcoxon test
p+stat_compare_means()

image

# 使用t-test统计检验方法
p+stat_compare_means(method = "t.test")

image

上述显著性标记可以通过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与显著性水平分行显示

举个栗子：

# 显示显著性水平，位置在1.5两组间和Y轴40位置
p+stat_compare_means(aes(label=..p.signif..), label.x = 1.5, label.y = 40)
# 也可以将标签指定为字符向量，不要映射，只需将p.signif两端的..去掉即可
p+stat_compare_means(label = "p.signif", label.x = 1.5, label.y = 40)

比较两个paired sample

# 比较两个paired sample
compare_means(len~supp, data=ToothGrowth, paired = TRUE)# 利用ggpaired()进行可视化
ggpaired(ToothGrowth, x="supp", y="len", color = "supp", line.color = "gray",line.size = 0.4, palette = "jco")+ stat_compare_means(paired = TRUE)

image

多组比较 Global test

anova进行多组比较

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.532727e-16 9.532727e-16  9.5e-16     ****  Anova

可视化, default Kruskal-Wallis

ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+stat_compare_means()

image

使用其他的方法，anova

#使用其他的方法, anova
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+stat_compare_means(method = "anova")

image

成对比较

默认进行两两比较

# 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   method<chr>  <chr>  <chr>        <dbl>        <dbl>    <chr>    <chr>    <chr>
1   len    0.5      1 7.020855e-06 1.404171e-05  7.0e-06     **** Wilcoxon
2   len    0.5      2 8.406447e-08 2.521934e-07  8.4e-08     **** Wilcoxon
3   len      1      2 1.772382e-04 1.772382e-04  0.00018      *** Wilcoxon

可以指定比较哪些组

# 可以指定比较哪些组
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

image

可以通过修改参数label.y来更改标签的位置

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

至于通过添加线条来连接比较的两组，这一功能已由包ggsignif实现

设定参考组

设定参考组

compare_means(len~dose, data=ToothGrowth, ref.group = "0.5",  #以dose=0.5组为参考组method = "t.test" )
# 可视化
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+stat_compare_means(method = "anova", label.y = 40)+ # Add global p-valuestat_compare_means(label = "p.signif", method = "t.test", ref.group = "0.5") # Pairwise comparison against reference

image

参考组也可以设置为.all.即所有的平均值

# 参考组也可以设置为.all.即所有的平均值
compare_means(len~dose, data=ToothGrowth, ref.group = ".all.", method = "t.test")
#可视化
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+stat_compare_means(method = "anova", label.y = 40)+# Add global p-valuestat_compare_means(label = "p.signif", method = "t.test",ref.group = ".all.")#Pairwise comparison against all

image

为什么有时需要将ref.group设置为.all

接下来利用survminer包中的数据集myeloma来讲解一下为什么有时候我们需要将ref.group设置为.all.

# 利用survminer包中的数据集myeloma来讲解一下为什么有时候我们需要将ref.group设置为.all.
install.packages("survminer")
library(survminer) #没安装的先安装再加载
data("myeloma")
head(myeloma)

我们将根据患者的分组来绘制DEPDC1基因的表达谱，看不同组之间是否存在显著性的差异，我们可以在7组之间进行比较，但是这样的话组间比较的组合就太多了，因此我们可以将7组中每一组与全部平均值进行比较，看看DEPDC1基因在不同的组中是否过表达还是低表达。

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<chr>  <chr>            <chr>        <dbl>        <dbl>    <chr>    <chr>
1 DEPDC1  .all.       Cyclin D-1 2.877529e-01 1.000000e+00     0.29       ns
2 DEPDC1  .all.       Cyclin D-2 4.244240e-01 1.000000e+00     0.42       ns
3 DEPDC1  .all.     Hyperdiploid 2.725486e-08 1.907840e-07  2.7e-08     ****
4 DEPDC1  .all. Low bone disease 5.258400e-06 3.155040e-05  5.3e-06     ****
5 DEPDC1  .all.              MAF 2.538126e-01 1.000000e+00     0.25       ns
6 DEPDC1  .all.            MMSET 5.784193e-01 1.000000e+00     0.58       ns
7 DEPDC1  .all.    Proliferation 2.393921e-05 1.196961e-04  2.4e-05     ****
# ... with 1 more variables: method <chr>

可视化DEPDC1基因表达谱

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 meanstat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-valuestat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.")# Pairwise comparison against all

image

从图中可以看出，DEPDC1基因在Proliferation组中显著性地过表达，而在Hyperdiploid和Low bone disease显著性地低表达

我们也可以将非显著性标记ns去掉，只需要将参数hide.ns=TRUE

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 meanstat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-valuestat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.", hide.ns = TRUE)# Pairwise comparison against all

分面多组同时比较

按另一个变量进行分组之后进行统计检验，比如按变量dose进行分组：

compare_means(len~supp, data=ToothGrowth, group.by = "dose")
#可视化
p <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp",palette = "jco", add = "jitter", facet.by = "dose", short.panel.labs = FALSE)#按dose进行分面
#label只绘制p-value
p+stat_compare_means(label = "p.format")

image

# label绘制显著性水平
p+stat_compare_means(label = "p.signif", label.x = 1.5)

image

将所有箱线图绘制在一个panel中

# 将所有箱线图绘制在一个panel中
p <- ggboxplot(ToothGrowth, x="dose", y="len", color = "supp",palette = "jco", add = "jitter")
p+stat_compare_means(aes(group=supp))
# 只显示p-value
p+stat_compare_means(aes(group=supp), label = "p.format")
# 显示显著性水平
p+stat_compare_means(aes(group=supp), label = "p.signif")

image

进行paired sample检验

compare_means(len~supp, data=ToothGrowth, group.by = "dose", paired = TRUE)
# 可视化
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) # 按dose分面
# 只显示p-value
p+stat_compare_means(label = "p.format", paired = TRUE)

image

其他图形

有误差棒的条形图

实际上我以前的文章里有纯粹用ggplot2实现

ggbarplot(ToothGrowth, x="dose", y="len", add = "mean_se")+stat_compare_means()+stat_compare_means(ref.group = "0.5", label = "p.signif", label.y = c(22, 29))

image

有误差棒的线图

ggline(ToothGrowth, x="dose", y="len", add = "mean_se")+stat_compare_means()+stat_compare_means(ref.group = "0.5", label = "p.signif", label.y = c(22, 29))

image

条形图(两个分组变量)

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)

image

线图-两分组

# line, multiply pair group
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))

image

环境信息

sessionInfo()

信息如下：

R version 3.4.3 (2017-11-30)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)Matrix products: defaultlocale:
[1] LC_COLLATE=Chinese (Simplified)_China.936 
[2] LC_CTYPE=Chinese (Simplified)_China.936   
[3] LC_MONETARY=Chinese (Simplified)_China.936
[4] LC_NUMERIC=C                              
[5] LC_TIME=Chinese (Simplified)_China.936    attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     other attached packages:
[1] survminer_0.4.2  bindrcpp_0.2     ggpubr_0.1.6.999 magrittr_1.5    
[5] ggplot2_2.2.1   loaded via a namespace (and not attached):[1] Rcpp_0.12.13        compiler_3.4.3      plyr_1.8.4         [4] bindr_0.1           tools_3.4.3         digest_0.6.15      [7] tibble_1.3.4        gtable_0.2.0        nlme_3.1-131       
[10] lattice_0.20-35     pkgconfig_2.0.1     rlang_0.1.4        
[13] Matrix_1.2-12       psych_1.7.8         ggsci_2.8          
[16] cmprsk_2.2-7        yaml_2.1.15         parallel_3.4.3     
[19] gridExtra_2.3       knitr_1.19          dplyr_0.7.4        
[22] stringr_1.2.0       survMisc_0.5.4      grid_3.4.3         
[25] tidyselect_0.2.3    data.table_1.10.4-3 glue_1.2.0         
[28] KMsurv_0.1-5        R6_2.2.2            km.ci_0.5-2        
[31] survival_2.41-3     foreign_0.8-69      tidyr_0.8.0        
[34] purrr_0.2.4         reshape2_1.4.2      splines_3.4.3      
[37] scales_0.5.0        assertthat_0.2.0    mnormt_1.5-5       
[40] xtable_1.8-2        colorspace_1.3-2    ggsignif_0.4.0     
[43] labeling_0.3        stringi_1.1.5       lazyeval_0.2.1     
[46] munsell_0.4.3       broom_0.4.3         zoo_1.8-1

猜你喜欢

  • 热文：1高分文章 2不可或缺的人 3图表规范
  • 一文读懂：1微生物组 2寄生虫益处 3进化树
  • 必备技能：1提问 2搜索 3Endnote
  • 文献阅读 1热心肠 2SemanticScholar 3geenmedical
  • 扩增子分析：1图表解读 2分析流程 3统计绘图 4功能预测
  • 科研经验：1云笔记 2云协作 3公众号
  • 系列教程：1Biostar 2微生物组 3宏基因组
  • 生物科普 1肠道细菌 2人体上的生命 3生命大跃进 4细胞的暗战 5人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1200+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
image

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
image

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

这篇关于R语言添加p-value和显著性标记的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/665700。

相关文章

Java枚举类实现Key-Value映射的多种实现方式

Java枚举类实现Key-Value映射的多种实现方式

《Java枚举类实现Key-Value映射的多种实现方式》在Java开发中,枚举(Enum)是一种特殊的类,本文将详细介绍Java枚举类实现key-value映射的多种方式,有需要的小伙伴可以根据需要... 目录前言一、基础实现方式1.1 为枚举添加属性和构造方法二、http://www.cppcns.co
阅读更多...
如何解决mysql出现Incorrect string value for column ‘表项‘ at row 1错误问题

如何解决mysql出现Incorrect string value for column ‘表项‘ at row 1错误问题

《如何解决mysql出现Incorrectstringvalueforcolumn‘表项‘atrow1错误问题》：本文主要介绍如何解决mysql出现Incorrectstringv... 目录mysql出现Incorrect string value for column ‘表项‘ at row 1错误报错
阅读更多...
C语言中的数据类型强制转换

C语言中的数据类型强制转换

《C语言中的数据类型强制转换》：本文主要介绍C语言中的数据类型强制转换方式,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录C语言数据类型强制转换自动转换强制转换类型总结C语言数据类型强制转换强制类型转换：是通过类型转换运算来实现的，主要的数据类型转换分为自动转换
阅读更多...
利用Go语言开发文件操作工具轻松处理所有文件

利用Go语言开发文件操作工具轻松处理所有文件

《利用Go语言开发文件操作工具轻松处理所有文件》在后端开发中,文件操作是一个非常常见但又容易出错的场景,本文小编要向大家介绍一个强大的Go语言文件操作工具库,它能帮你轻松处理各种文件操作场景... 目录为什么需要这个工具？核心功能详解1. 文件/目录存javascript在性检查2. 批量创建目录3. 文件
阅读更多...
C语言实现两个变量值交换的三种方式

C语言实现两个变量值交换的三种方式

《C语言实现两个变量值交换的三种方式》两个变量值的交换是编程中最常见的问题之一,以下将介绍三种变量的交换方式,其中第一种方式是最常用也是最实用的,后两种方式一般只在特殊限制下使用,需要的朋友可以参考下... 目录1.使用临时变量（推荐）2.相加和相减的方式（值较大时可能丢失数据）3.按位异或运算1.使用临时
阅读更多...
使用C语言实现交换整数的奇数位和偶数位

使用C语言实现交换整数的奇数位和偶数位

《使用C语言实现交换整数的奇数位和偶数位》在C语言中,要交换一个整数的二进制位中的奇数位和偶数位,重点需要理解位操作,当我们谈论二进制位的奇数位和偶数位时,我们是指从右到左数的位置,本文给大家介绍了使... 目录一、问题描述二、解决思路三、函数实现四、宏实现五、总结一、问题描述使用C语言代码实现：将一个整
阅读更多...
C语言字符函数和字符串函数示例详解

C语言字符函数和字符串函数示例详解

《C语言字符函数和字符串函数示例详解》本文详细介绍了C语言中字符分类函数、字符转换函数及字符串操作函数的使用方法,并通过示例代码展示了如何实现这些功能,通过这些内容,读者可以深入理解并掌握C语言中的字... 目录一、字符分类函数二、字符转换函数三、strlen的使用和模拟实现3.1strlen函数3.2st
阅读更多...
Go语言中最便捷的http请求包resty的使用详解

Go语言中最便捷的http请求包resty的使用详解

《Go语言中最便捷的http请求包resty的使用详解》go语言虽然自身就有net/http包,但是说实话用起来没那么好用,resty包是go语言中一个非常受欢迎的http请求处理包,下面我们一起来学... 目录安装一、一个简单的get二、带查询参数三、设置请求头、body四、设置表单数据五、处理响应六、超
阅读更多...
C语言中的浮点数存储详解

C语言中的浮点数存储详解

《C语言中的浮点数存储详解》：本文主要介绍C语言中的浮点数存储详解,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录1、首先明确一个概念2、接下来，讲解C语言中浮点型数存储的规则2.1、可以将上述公式分为两部分来看2.2、问：十进制小数0.5该如何存储？2.3 浮点
阅读更多...
基于Python实现多语言朗读与单词选择测验

基于Python实现多语言朗读与单词选择测验

《基于Python实现多语言朗读与单词选择测验》在数字化教育日益普及的今天,开发一款能够支持多语言朗读和单词选择测验的程序,对于语言学习者来说无疑是一个巨大的福音,下面我们就来用Python实现一个这... 目录一、项目概述二、环境准备三、实现朗读功能四、实现单词选择测验五、创建图形用户界面六、运行程序七、
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2