厚缊 业余的R语言可视化重度患者 个人博客:houyun.xyz 邮箱:houyunhuang@163.com 转载本文(包括长期转载账号)必须联系厚缊授权 相关系数矩阵可视化已经至少有两个版本的实现了,魏太云基于base绘图系统写了corrplot 包,应该说是相关这个小领域中最精美的包了,使用简单,样式丰富,只能用惊艳来形容。Kassambara的ggcorrplot 基于ggplot2 重写了corrplot ,实现了corrplot 中绝大多数的功能,但仅支持“square”和“circle”的绘图标记,样式有些单调,不过整个ggcorrplot 包的代码大概300行,想学习用ggplot2 来自定义绘图函数,看这个包的源代码很不错。还有部分功能相似的corrr 包(在写ggcor 之前完全没有看过这个包,写完之后发现在相关系数矩阵变data.frame方面惊人的相似),这个包主要在数据相关系数提取、转换上做了很多的工作,在可视化上稍显不足。ggcor 的核心是为相关性分析、数据提取、转换、可视化提供一整套解决方案,目前的功能大概完成了70%,后续会根据实际需要继续扩展。这里要介绍的ggcor 是corrplot 的有一种实现,在吸收借鉴(或者说是全般)corrplot 的基础上,略有提升,使用上会更灵活简单。这里不得不说下整个开发过程,一方面是为了感谢来自五湖四海的朋友提的许许多多的专业意见,另一方面是深感歉意,由于我考虑不周,ggcor (第一次上传github测试大约在十天前,名字还是ggcorrr )存在各种各样的小bug,也有很多限制,这次为了升级和简化,不仅名字变了,里面的函数也变得完全不一样,若你们以前写了ggcorrr 的代码,现在基本不可以重用了。尽管这一版本的更新我考虑了很多应用场景,综合了方方面面的意见,但是bug和限制还是在所难免,希望朋友们继续提意见,提需求。当然,这一次有两点是可以保证的:一是包名不会变了,再变都成精了;二是包内的主体函数和参数不会变了,而且即使要变,我也会考虑兼容性。安装if(!require(devtools)) install.packages('devtools') if(!require(ggcor)) devtools::install_github('houyunhuang/ggcor')
数据预处理函数ggplot2 要求的数据格式是data.frame ,要把相关系数矩阵处理成理想的数据格式需要一系列的操作。ggcor 提供了两个主要的函数,一个是as_cor_tbl() 函数,另一个是fortify_cor() 函数,两者结合应该能满足绝大多数的应用场景需求。as_cor_tbl() 是更底层的函数,fortify_cor() 本质上调用了as_cor_tbl() 来得到最终的数据结果。这个函数适用于已经知道,或者需要用其它更特殊的函数(非stats::cor() )来处理得到系数的情况,常用的参数是前三个。x —— 相关系数矩阵(或者数据框),矩阵行名和列名是必要的,若没有或者缺失值会自动补全名字,行名以“Y”开头,附上递增的整数序列,列名以“X”开头,附上附上递增的整数序列。type —— 相关系数矩阵图样式,“upper”截断下三角,“lower”截断上三角。show.diag —— 相关系数矩阵图中是否包含对角线,仅对对称矩阵有效。p —— 相关系数检验p值矩阵(或者数据框),必须与x 一一对应。low —— 相关系数置信区间下界矩阵(或者数据框),必须与x 一一对应。upp —— 相关系数置信区间上界矩阵(或者数据框),必须与x 一一对应。cluster.type —— 是否对相关系数矩阵进行重新排序,“none”表示不重排,“all”表示行列均重排,“row”表示对行重排,“col”则只对列重排。... —— 其它传递给matrix_order() 函数的参数。
library(ggcor) ## function(x, ## type = c('full', 'upper', 'lower'), ## show.diag = TRUE, ## p = NULL, ## low = NULL, ## upp = NULL, ## cluster.type = c('none', 'all', 'row', 'col'), ## ...)
corr <- cor(mtcars) df <- as_cor_tbl(corr) df ## return a tibble
## # A tibble: 121 x 3 ## x y r ## * <int> <int> <dbl> ## 1 1 11 1 ## 2 1 10 -0.852 ## 3 1 9 -0.848 ## 4 1 8 -0.776 ## 5 1 7 0.681 ## # … with 116 more rows ## Extra attributes: ## xname: mpg cyl disp hp drat wt qsec vs am gear carb ## yname: carb gear am vs qsec wt drat hp disp cyl mpg ## show.diag: TRUE
fortify_cor() 主要适用于处理原始数据表,即调用cor() 求相关系数,cor 函数对数据按列进行两两相关性计算,默认使用pearson 方法,当然理论解读中提前的spearman 和kendall 方法也都支持。调用(若需要)cor.test() 函数进行统计检验,并返回ggcor 需要的数据类。x —— 原数据矩阵(或者数据框),列名是必要的,若没有或者缺失值会自动补全名字,列名以“X”开头,附上附上递增的整数序列。y —— 原数据矩阵(或者数据框),列名是必要的,若没有或者缺失值会自动补全名字,列名以“X”开头,附上附上递增的整数序列。当y 不为空(NULL)时,相关系数是x 中的每一列和y 中的每一列的相关性。type —— 相关系数矩阵图样式,“upper”截断下三角,“lower”截断上三角。show.diag —— 相关系数矩阵图中是否包含对角线,仅对对称矩阵有效。cor.test —— 逻辑值,是否进行相关性检验。cor.test.alt —— 相关性检验备择假设,详细请查看cor.test() 帮助。cor.test.method —— 相关性检验方法,详细请查看cor.test() 帮助。cluster.type —— 是否对相关系数矩阵进行重新排序,“none”表示不重排,“all”表示行列均重排,“row”表示对行重排,“col”则只对列重排。cluster.method —— 当cluster.order 为“HC”(默认)时算法,详细请查看ggcor::matrix_order() 。
## function( ## x, ## y = NULL, ## type = c('full', 'upper', 'lower'), ## show.diag = FALSE, ## cor.test = FALSE, ## cor.test.alt = 'two.sided', ## cor.test.method = 'pearson', ## cluster.type = c('none', 'all', 'row', 'col'), ## cluster.method = 'HC', ## ... )
df01 <- fortify_cor(mtcars, cor.test = TRUE, cluster.type = 'all') df01
## # A tibble: 121 x 6 ## x y r p low upp ## * <int> <int> <dbl> <dbl> <dbl> <dbl> ## 1 1 11 1 7.44e-232 1.000 1.000 ## 2 1 10 0.750 7.83e- 7 0.543 0.871 ## 3 1 9 0.428 1.46e- 2 0.0927 0.676 ## 4 1 8 0.527 1.94e- 3 0.218 0.740 ## 5 1 7 0.395 2.53e- 2 0.0537 0.654 ## # … with 116 more rows ## Extra attributes: ## xname: carb hp wt cyl disp drat am gear qsec mpg vs ## yname: vs mpg qsec gear am drat disp cyl wt hp carb ## show.diag: TRUE
as_cor_tbl() 和fortify_cor() 返回值均是cor_tbl 类的数据框(准确的说是tibble ),并包含其它额外特殊属性,要充分利用ggcor 包中一系列工具函数带来的便捷性,必须调用(手动比较麻烦)这两个函数预处理数据。初始化绘图函数ggcor 封装了一个基本的初始化函数ggcor() ,用来处理数据、绘图类型、背景、坐标轴、颜色映射、图例等。看上去ggcor() 函数非常复杂,但若仔细观察参数命名规则和分类,又是非常简单的。绘图区网格线参数名字均以grid.* 开头;坐标轴相关参数均以axis.* 开头;图例相关(主要是colourbar)的参数均以legend.* 开头;还有panel.backgroud 、xlim 、ylim 均是常见的参数,panel.backgroud 参数用来设置绘图区背景颜色,xlim 、ylim 则是设置x/y轴的范围。把这几大块参数去掉后,剩下的参数只有6个了,瞬间清爽了不少。要完全理解ggcor() 的作用原理及相关参数的设置,需要先讲讲ggcor() 内部构成。ggcor() 本质上是调用了ggplot() 来初始化,然后根据相关系数图的样式添加了一些辅助的图层。x 、y 、mapping 、is.cor 、show.diag 和... 参数均和数据预处理和映射相关。x 可以是cor_tbl 、矩阵(数据框)。当为cor_tbl 时直接作为data 参数传递给ggplot() ;为矩阵(数据框)时,若是(is.cor = TRUE )相关系数矩阵(数据框)时,调用as_cor_tbl() 函数处理成cor_tbl ,若不是(is.cor = FALSE )相关系数矩阵(数据框)时,调用fortify_cor() 函数处理成cor_tbl ,此时x 、y 、show.diag 和... 均传递给as_cor_tbl() 或者fortify_cor() 。mapping 对应ggplot() 中的mapping 参数,当为空(默认)时,根据cor_tbl 中的变量情况添加,基础形式是aes(x = x, y = y, r = r, fill = r) 。若cor_tbl 包含“p”(进行了相关系数显著性检验),则最基础形式基础上额外添加p = p ,若检验方法(cor.test.method = 'pearson' ),再加上low = low 和upp = upp 。
- 若
fill.scale.add 为FALSE 不添加颜色映射函数。若为TRUE (默认),则会在初始化中自动添加颜色映射函数。 fill.colours 是颜色字符向量,用来控制填充色(fill)的颜色映射,若为空(默认),使用corrplot 默认配色,可以通过ggcor:::.defualt.colors 查看具体颜色。fill.bin 是否对连续性颜色进行分组(默认FALSE )。当fill.scale.add = TRUE 时,若fill.bin = TRUE ,ggcor() 的颜色映射函数是使用scale_fill_steps2n() ,若fill.bin = FALSE ,ggcor() 的颜色映射函数是使用scale_fill_gradient2n() 。
- 背景网格线是通过
geom_tile() (目前来看,geom_tile() 会存在一些难以处理的缺点,以后可能会添加专用的网格线图层函数),panel.backgroud (fill)、grid.* 系列参数传递给geom_tile() 。 legend.position 传递给theme() 中对应的参数,用来控制图例位置,其它legend.* 开头参数传递给guide_colourbar() 或者guide_colorsteps() 。legend.breaks 用来控制图例颜色棒标签显示位置,legend.labels 是对应的标签。若fill.bin = TRUE ,legend.breaks 也是图例颜色棒切割分组的位置。coord.fixed 逻辑值,为真xlim 、ylim 传递给coord.fixed() 函数,为假传递给coord_cartesian() 函数。在ggcor 包中,相关系数矩阵若是n * m的矩阵,那么第i行对应的坐标点(即as_cor_tbl() 返回结果中的y)为n-i(为了和表格呈现样式一致,行方向翻转了),第j列对应的坐标点(即as_cor_tbl() 返回结果中的x)为j,得到第(i, j)个数据点所在的方格坐标为(xmin = j-0.5, xmax = j+0.5, ymin = n-i-0.5, ymax = n-i+0.5)。从这个规律我们不难算出默认的xlim 、ylim 取值范围。注意,对cor_tbl 数据框,或者说是ggcor 包中函数操作数据,均不会改变每个数据单元格在图中的坐标位置。
ggcor() 初始化之后,本质上返回的是ggplot 对象,若是想改变默认设置,可以按照ggplot2 的相应的函数和设置方法去调整。## function( ## x, ## y = NULL, ## mapping = NULL, ## is.cor = FALSE, ## show.diag = TRUE, ## fill.scale.add = TRUE, ## fill.colours = NULL, ## fill.bin = FALSE, ## panel.backgroud = NA, ## grid.colour = 'grey50', ## grid.size = 0.25, ## grid.linetype = 'solid', ## axis.x.position = c('auto', 'bottom', 'top'), ## axis.y.position = c('auto', 'left', 'right'), ## fill.colours = NULL, ## legend.title = 'correalation', ## legend.position = 'auto', ## legend.breaks = NULL, ## legend.labels = NULL, ## coor.fixed = TRUE, ## xlim = NULL, ## ylim = NULL, ## ...) 我写ggcor() 函数花的时间最长,耗费的精力最多,但到现在仍然是半成品的感觉,有些问题直到现在我也没找到相对完美的解决方案。尽管如此,对于新手,我还是建议调用ggcor() 来进行初始化,若自己去研究各种图层函数,折腾很多细节,一天也难得出一幅图,对于心里的打击比较大。看几个初始化之后的效果。ggcor(mtcars)
ggcor(mtcars, type = 'lower')
df02 <- fortify_cor(mtcars, type = 'upper') ggcor(df02, panel.backgroud = '#66C2A5')
图层函数ggcor 提供了定制的geom_square() 、geom_circle2() 、geom_ellipse2() 、geom_pie2() 、geom_colour() 、geom_confbox() 、geom_num() 、geom_mark() 、geom_cross() 9个ggplot2 图层函数,可以根据需要进行叠加。除了ggplot2 中一般化的参数(x 、y 、fill 、colour 、size 等)最常用参数r 、p 、low 、upp 、num 、r0 、sig.thres 、sig.level 、mark 等。r —— 相关系数,适用于geom_square() 、geom_circle2() 、geom_ellipse2() 、geom_pie2() 、geom_confbox() 、geom_mark() (统计显著性标记的系数,与 1.234∗∗类似)。这些参数之所以都设置为“r”,主要是因为在相关系数可视化中基本都映射为相关系数,统一命名可以减少一些参数记忆,方便使用。p —— 相关系数检验P值,适用于geom_mark() 、geom_cross() ,结合sig.thres 等参数来根据显著性水平做一些辅助标记。r 、low 、upp —— 适用于geom_confbox() ,三个参数均是必须的,low 、r 、upp 分别确定置信区间盒子的下端、中间和上端线条位置。num —— 适用于geom_num() ,数值变量。这个函数封装于geom_text() ,做了一点点优化,以后可能会删除。
r0 ,外接圆半径缩放系数,适用于geom_square() 、geom_circle2() 、geom_ellipse2() 、geom_pie2() 函数。该参数的主要意义是处理图形覆盖问题,当在每个单元格画半径为0.5的方块、圆等图标时,会相互覆盖掉背景网格线,影响视觉效果。该参数默认值是0.48。sig.thres 统计显著性临界值,用来过滤非显著的数据行。适用于geom_mark() 、geom_cross() 。sig.level 、mark 适用于geom_mark() ,前者为统计显著性水平向量(如c(0.001, 0.01, 0.05) ),后者为对应的标记符号(c('***', '**', '*') )。统计显著性水平向量不要求一定要按顺序排列,只要求和标记符号一一对应就行。
ggcor(mtcars) + geom_square()
ggcor(mtcars, type = 'upper') + geom_circle2()
ggcor(mtcars, type = 'lower', show.diag = TRUE) + geom_ellipse2()
ggcor(mtcars, type = 'full', cluster.type = 'all') + geom_pie2()
ggcor(mtcars, cluster.type = 'all') + geom_colour() + geom_num(aes(num = r), colour = 'grey90', size = 3.5)
ggcor(mtcars, type = 'full', cor.test = TRUE) + geom_confbox()
ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_colour() + geom_cross()
ggcor(mtcars, type = 'full', cor.test = TRUE) + geom_square() + geom_cross() 其它图层函数的使用都比较符合直觉,需要设置的地方也很少,而geom_mark() 会涉及一些其它问题,这里拉出来说说。很多时候,我们并不关心不具备统计显著性的相关系数,也不需要在图中显示,这时需要设置sig.thres ,即要过滤的显著性临界值。ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_raster() + geom_mark(sig.thres = 0.05, size = 3, colour = 'grey90') 若是只要统计显著性标记的'*'号,不要系数怎么整?在geom_mark() 中也很简单,直接设置r 就好。ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_raster() + geom_mark(r = NA, sig.thres = 0.05, size = 5, colour = 'grey90')
那要改变星号标记规则,只要小于0.05的标记一颗星,其它什么都不标记呢?注意:因为星号在文本中显示在偏上的位置,若不设置vjust 参数,看上去纵向会不居中。ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_raster() + geom_mark(r = NA, sig.level = 0.05, mark = '*', vjust = 0.65, size = 6, colour = 'grey90')
非对称相关系数矩阵非对称相关系数矩阵和非对称矩阵是有细微的区别的,前者表示行列代表不同的变量集合,相互之间的顺序可以打乱。所以,有时候要分析两个表中每个变量之间的相关性,此时得到的结果就是非对称的相关系数矩阵。library(vegan) # 使用vegan包所带的数据集 data(varechem) data(varespec) df03 <- fortify_cor(x = varechem, y = varespec[ , 1:30], cluster.type = 'col') ggcor(df03) + geom_colour()
df04 <- fortify_cor(x = varespec[ , 1:30], y = varechem, cor.test = TRUE) ggcor(df04) + geom_square() + geom_cross(size = 0.2)# 数据集不好,没几个显著的
上下三角不一样怎么画?ggcor 提供了很多辅助函数来对cor_tbl 数据进行过滤,函数命名规则上都以get_* 开头。get_lower_data() —— 获取相关系数矩阵下三角所在行,仅支持对称的相关系数矩阵。get_upper_data() —— 获取相关系数矩阵上三角所在行,仅支持对称的相关系数矩阵。get_diag_data() —— 获取相关系数矩阵对角线所在行,仅支持对称的相关系数矩阵。get_diag_tri() —— 删除相关系数矩阵对角线所在行,仅支持对称的相关系数矩阵。get_data() —— 是以上四个函数的重新包装,主要在画图时使用,稍后通过例子详细说明。
df05 <- fortify_cor(x = varechem, cor.test = TRUE, cluster.type = 'all') ggcor(df05) + geom_circle2()
df05_lower <- get_lower_data(df05, show.diag = FALSE) ggcor(df05_lower) + geom_circle2()
ggcor(df05) + geom_pie2(data = get_data(type = 'upper', show.diag = FALSE)) + geom_ellipse2(data = get_data(type = 'lower', show.diag = TRUE))
ggcor(df05) + geom_segment(aes(x = x - 0.5, y = y + 0.5, xend = x + 0.5, yend = y - 0.5), data = get_data(type = 'diag'), size = 0.5, colour = 'grey60') + geom_colour(data = get_data(type = 'upper', show.diag = FALSE)) + geom_mark(data = get_data(type = 'upper', show.diag = FALSE), size = 3) + geom_circle2(data = get_data(r >= 0.5, type = 'lower', show.diag = FALSE), r = 0.8, fill = '#66C2A5') + geom_num(aes(num = r), data = get_data(type = 'lower', show.diag = FALSE), size = 3)
列名放在对角线?
ggcor(mtcars, cor.test = TRUE, cluster.type = 'all') + geom_confbox(data = get_data(type = 'upper', show.diag = FALSE)) + geom_num(aes(num = r), data = get_data(type = 'lower', show.diag = FALSE), size = 3.5) + add_diaglab(size = 4.56) + remove_axis()
想对颜色分组?
很多情况下,连续性颜色棒并不是很好分区每个单元格对应的数值区间,这时根据相关系数大小对颜色进行分组可能更适合。说来非常巧合,一直头疼这个问题难以解决,就在前不久,ggraph 包的作者Thomas Lin向ggplot2 包提交了一组新的颜色映射函数(scale_fill/colour_steps*() ),这个问题迎刃而解,非常幸运。ggcor() 函数有参数fill.binned ,默认为FALSE ,设置为TRUE 就会根据相关系数大小对颜色分组。若要控制分组的数量和区间,可以通过legend.breaks 来设置。ggcor(mtcars, fill.bin = TRUE) + geom_square() # 默认分组
ggcor(mtcars, fill.bin = TRUE, legend.breaks = seq(-1, 1, length.out = 11)) + geom_square() #指定分组,0.2为一个区间
col <- col <- ggcor:::.default_colors col[6] <- '#F2F2F2' ggcor(mtcars, cluster.type = 'all', fill.colours = col, fill.bin = T, legend.breaks = c(-1, -0.8, -0.5, 0.5, 0.8, 1)) + geom_colour()
玩点花活
这部分内容要在线下载表情,很多时候会因为网络问题下载失败。不给图了。library(ggimage) emoji <- c('1f004', '1f0cf', '1f170', '1f171', '1f17e', '1f17f', '1f18e', '1f191', '1f192', '1f193', '1f194', '1f195', '1f196', '1f197') ggcor(df05) + geom_pokemon(aes(image=ifelse(r > 0.5, 'pikachu', 'tauros')), data = get_data(type = 'lower', show.diag = FALSE)) + geom_emoji(aes(image = ifelse(p <= 0.05, '1f600', '1f622')), data = get_data(type = 'upper', show.diag = FALSE)) + geom_emoji(aes(image = emoji), data = get_data(type = 'diag'))
ggcor(df05) + geom_pokemon(aes(image=ifelse(r > 0.5, 'pikachu', 'tauros')), data = get_data(type = 'lower', show.diag = FALSE)) + geom_colour(data = get_data(type = 'upper', show.diag = FALSE)) + geom_shade(data = get_data(type = 'upper', show.diag = FALSE), sign = -1, size = 0.1) + geom_emoji(aes(image = emoji), data = get_data(type = 'diag'))
mantel 检验组合图
mantel 检验(Mantel test 是对两个矩阵相关关系的检验)的组合图已经十分流行了,用各种工具做的都有。大概5月份的时候,我基于corrplot 模拟重现了那幅图,直到现在每周都有人询问我相关实现的问题,我基本都是回答说等新方案,因为那个实现很复杂,没有基本的R知识,很难替换成自己的数据。8月底,我在ggcorrr (ggcor 前身)开发了gghpcc ,几乎可以满足可视化的要求,但是在优化ggcor 的过程中,我才恍然大悟,mantel检验本质上也是相关性分析,何不统一整合进ggcor 呢,就这样,gghpcc 被干掉了。进行案例展示之前,我先解释一个关键性问题,这个问题在目前框架下没有很好的办法去自动化解决方案,那就是坐标轴范围,需要手动设置。前文有提及,这里再重复一次,ggcor() 初始化的默认坐标范围是和correlation matrix的行列数相关的,若行数为n,列数为m,x轴范围是c(0.5, m + 0.5),y轴的范围是c(0.5, n + 0.5)。数据预处理函数ggcor 提供了mantel检验的封装函数fortify_mantel() ,支持vegan 包中的mantel() 、mantel.partial() 和ade4 包中的mantel.randtest() 、mantel.rtest() 函数,差别上说mantel.partial() 是偏mantel检验(有控制变量),其它三个是mantel检验,当不使用并行计算时,mantel.randtest() 速度最快(底层是C语言),mantel.rtest() 最慢,纯粹R代码实现。spec 是物种数据,支持列表(list,非data.frame)或者数据框(data.frame)。- 若是数据框(最常见的情况),数据框中的每一列是一个物种(OTU),每行是一个样本,可以通过
spec.select 参数来指定哪些列构成一个群落。例如若spec 中1-4列是spec01群落,5-12是spec02群落,spec.select = list(spec01 = 1:4, spec02 = 5:12) (当然,你也可以(最好不)不指定群落名称,这是名字自动命名)。还有一种情况(设置spec.group 参数的情况)后面单独解释。
env 是环境数据,支持列表(list,非data.frame)或者数据框(data.frame),env 中的每个元素对应一个环境变量(当然,若是列表,也可以支持多个环境变量组合成一个环境因素的情况)。还有一种情况(设置env.group 参数的情况)后面单独解释。env.ctrl 是环境控制变量,支持列表(list,非data.frame)或者数据框(data.frame)。需要注意,当env.ctrl 非列表时,每次计算的控制环境是相同的,若需要分别设置不同的控制环境,需要通过列表手动设置。还有一种情况(设置env.ctrl.group 参数的情况)后面单独解释。mantel.fun mantel函数,“mantel”、“mantel.randtest”、“mantel.rtest”和“mantel.partial”。spec.dist.fun 物种距离函数,vegdist 或者dist 。env.dist.fun 环境距离函数,vegdist 或者dist 。注意,当前的情况下,控制环境使用环境距离函数。spec.dist.method 物种距离计算方法,参数默认是“bray”。env.dist.method 环境距离计算方法,参数默认是“euclidean”。 group相关的参数是为了处理需要根据样本进行分组的情况,比如我A、B、C三个不同的样本分组,物种、环境和控制环境(均必须为数据框)同样如此,可以通过向量索引(和样本量等长)来指定分组。这时,每个样本组的物种只和对应样本组的环境列表的每个元素进行mantel测试。注意:以上情况均需要设置is.pair = TRUE 。## function(spec, ## env, ## env.ctrl = NULL, ## mantel.fun = 'mantel', ## is.pair = FALSE, ## spec.select = NULL, # a list of index vector ## spec.group = NULL, ## env.group = NULL, ## env.ctrl.group = NULL, ## spec.dist.fun = 'vegdist', ## env.dist.fun = 'vegdist', ## spec.dist.method = 'bray', ## env.dist.method = 'euclidean', ## ...)
案例一:按列设置物种群落library(vegan) # 使用vegan包所带的数据集 data(varechem) data(varespec) fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20))
## # A tibble: 56 x 4 ## spec env r p ## <chr> <chr> <dbl> <dbl> ## 1 spec01 N 0.100 0.129 ## 2 spec02 N 0.254 0.011 ## 3 spec03 N 0.0596 0.238 ## 4 spec04 N 0.182 0.038 ## 5 spec01 P 0.165 0.028 ## 6 spec02 P 0.112 0.129 ## 7 spec03 P -0.0240 0.573 ## 8 spec04 P 0.201 0.02 ## 9 spec01 K 0.0481 0.295 ## 10 spec02 K 0.307 0.008 ## # … with 46 more rows
案例二:按样本分组group <- rep(LETTERS[1:3], 8) fortify_mantel(varespec[ , 38:43], varechem, spec.group = group, env.group = group, is.pair = TRUE, mantel.fun = 'mantel.randtest')
## # A tibble: 42 x 4 ## spec env r p ## <chr> <chr> <dbl> <dbl> ## 1 A N -0.334 0.884 ## 2 B N 0.212 0.167 ## 3 C N 0.241 0.223 ## 4 A P -0.295 0.888 ## 5 B P 0.250 0.099 ## 6 C P 0.312 0.14 ## 7 A K 0.0909 0.31 ## 8 B K 0.310 0.101 ## 9 C K 0.0720 0.33 ## 10 A Ca -0.365 0.942 ## # … with 32 more rows
案例三:偏mantel测试fortify_mantel(varespec, varechem[ , 1:10], env.ctrl = varechem[ , 11:14], spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest')
## # A tibble: 40 x 4 ## spec env r p ## <chr> <chr> <dbl> <dbl> ## 1 spec01 N 0.100 0.134 ## 2 spec02 N 0.254 0.017 ## 3 spec03 N 0.0596 0.243 ## 4 spec04 N 0.182 0.039 ## 5 spec01 P 0.165 0.031 ## 6 spec02 P 0.112 0.122 ## 7 spec03 P -0.0240 0.582 ## 8 spec04 P 0.201 0.009 ## 9 spec01 K 0.0481 0.271 ## 10 spec02 K 0.307 0.007 ## # … with 30 more rows
mantel检验可视化mantel检验既然是相关性分析的一种特殊情况,那么相关性分析的可视化也同样适用于mantel检验。当然,由于fortify_mantel() 返回数据类和fortify_cor() 返回的数据类并不一样,需要经过一道转换程序。mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20)) df06 <- as_cor_tbl(mantel) ggcor(df06) + geom_pie2() + geom_cross()
mantel检验组合图
mantel组合图是与相关性分析高度整合的,依赖于相关性分析函数,换句话说mantel组合图只是在相关性分析图的基础上额外叠加了一个图层。核心函数是add_link() ,是个泛型函数(后面另说)。参数暂时不解释了,看例子。实在抱歉,add_mantel() 又被我干掉了。corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest') ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
corr <- fortify_cor(varechem, type = 'upper', show.diag = FALSE, cor.test = TRUE, cluster.type = 'all') ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_pie2() + remove_axis('y')
corr <- fortify_cor(varechem, type = 'lower', show.diag = FALSE, cor.test = TRUE, cluster.type = 'all') ggcor(corr, xlim = c(0.5, 20)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_ellipse2() + remove_axis('y')
## 仅是测试效果,没有实际意义 corr <- fortify_cor(varechem, varechem[ , 1:7], type = 'full', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest', nrepet = 2000) extra.params <- extra_params(group.label = text_params(size = 6), link.params = link_params(group.point.hjust = 2)) ggcor(corr, axis.y.position = 'left', legend.position = 'left', xlim = c(0.5, 14.5)) + add_link(mantel, extra.params = extra.params) + geom_circle2()
group <- rep(LETTERS[1:3], 8)
corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec[ , 38:43], varechem, spec.group = group, env.group = group, is.pair = TRUE, mantel.fun = 'mantel.randtest')
ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_colour() + geom_shade(sign = -1) + remove_axis('y')
corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest') mantel <- dplyr::filter(mantel, p <= 0.05)
ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE, legend.drop = TRUE) + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
## 个人觉得很丑
只要简单的相关性?corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') corr01 <- fortify_cor(varechem, varespec[ , 38:39], type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest')
ggcor(corr, xlim = c(-5, 14.5)) + add_link(x = corr01, diag.label = TRUE, link.line.colours = c('#E31A1C', '#33A02C')) + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
ggcor(corr, xlim = c(-5, 14.5)) + add_link(x = corr01, mapping = aes(size = abs(r)), diag.label = TRUE) + add_diaglab(angle = 45) + geom_square() + scale_size_continuous(limits = c(0, 1), range = c(0.25, 3)) + guides(size = guide_legend(title = 'abs r', override.aes = list(colour = 'grey35'), order = 1)) + remove_axis('y')
更灵活的方案若要自己更灵活的处理,或者是要用这个组合图展示其它类型的数据,add_link() 支持自定义的数据框做参数。第一个参数df 需要一个数据框,包含x和y列,x列类似于mantel检验中的物种群落(或者是样本组),y类似于mantel检验中的环境变量。在内部,会自动用df 中的y和相关系数矩阵的行名进行匹配确定坐标位置。corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') df <- data.frame(x = rep(LETTERS[1:3], 14), y = rep(cor_tbl_yname(corr), 3), r = runif(42, -1, 1), p = runif(42, 0, 0.5), stringsAsFactors = FALSE) ggcor(corr, xlim = c(-5, 14.5)) + add_link(df, diag.label = TRUE, colour = 'red') + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
相关性网络图这块内容不会整合在ggcor 包里面,但是利用ggcor 里面的函数很容易导出相关性分析数据供其它函数使用。df07 <- fortify_cor(varespec, type = 'upper', show.diag = FALSE, cor.test = TRUE, keep.name = TRUE) df07
## # A tibble: 946 x 6 ## x y r p low upp ## <fct> <fct> <dbl> <dbl> <dbl> <dbl> ## 1 Empenigr Callvulg -0.262 0.217 -0.602 0.159 ## 2 Rhodtome Callvulg -0.0608 0.778 -0.453 0.351 ## 3 Rhodtome Empenigr 0.529 0.00781 0.160 0.769 ## 4 Vaccmyrt Callvulg -0.0925 0.667 -0.478 0.323 ## 5 Vaccmyrt Empenigr 0.205 0.337 -0.216 0.562 ## # … with 941 more rows ## Extra attributes: ## xname: Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv ... ## yname: Cladphyl Claddefo Cladcerv Icmaeric Peltapht Stersp N ... ## show.diag: FALSE
## make graph library(ggraph) library(tidygraph) graph <- as_tbl_graph(df07) ggraph(graph, layout = 'linear', circular = TRUE) + geom_edge_arc(aes(colour = r, alpha = p)) + scale_edge_alpha_continuous(range = c(1, 0.1)) + coord_fixed() #颜色是相关性,线条浓淡是统计检验P值
写在最后的话
后续除了优化现有的功能,最重要的一个模块是写一个交互界面,能让不懂R的人,或者不喜欢敲代码的人也能轻松出图。当然,任何好用的包,都是无数的人共同协作、共同努力的结果。我希望能有更多的人参与测试,更多的人提供一些与相关性分析相关的新功能。ggcor 还有很多的小技巧没有来得及介绍,全部留在下一期吧。高颜值免费在线绘图R统计和作图
|