厚缊

诹图——ggcor简介(一)

厚缊 / 2020-01-31


ggcor前前后后经过了好几轮更新,文档一直还是最初的样子,这次趁着疫情紧张,全国人民都在闭关做贡献的时间里,集中更新介绍文章。这次我会尽可能说得详细一些,可能会用5-7篇短文章来进行介绍。

安装

## install.packages("devtools")
devtools::install_github("houyunhuang/ggcor") 
## 若原来安装过老版本
devtools::install_github("houyunhuang/ggcor", force = TRUE) 

本系列文章都是基于ggcor0.9.2.11版本,若不是最新的请更新包。

packageVersion("ggcor")
## [1] '0.9.2.11'

相关系数计算和检验

既然是相关性分析和可视化,最核心的当然是相关系数的计算和检验,R 有两个最基本的函数:cor()计算相关系数矩阵,cor.test()进行相关系数统计检验。遗憾的是,cor.test()每次只能检验两个变量,并不能直接返回相关系数检验矩阵,所以需要做一定的封装来完成此项任务,这就是ggcor自带correlate()函数的用途。

correlate()函数

correlate()ggcor包中相关系数计算和检验的默认函数,了解这个函数有助于理解默认情况下ggcor是如何计算的。

library(ggcor)
args(correlate)
## function (x, y = NULL, cor.test = FALSE, method = "pearson", 
##     use = "everything", ...) 
## NULL

默认情况下,correlate()函数不会进行相关系数检验的。

correlate(mtcars) %>% str()
## List of 4
##  $ r       : num [1:11, 1:11] 1 -0.852 -0.848 -0.776 0.681 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ p.value : NULL
##  $ lower.ci: NULL
##  $ upper.ci: NULL
##  - attr(*, "class")= chr "correlation"

cor.test参数设置为TRUE,进行统计检验。method = "pearson"时,会返回统计检验置信区间矩阵,即返回值中的upper.ci和lower.ci。

correlate(mtcars, cor.test = TRUE) %>% str()
## List of 4
##  $ r       : num [1:11, 1:11] 1 -0.852 -0.848 -0.776 0.681 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ p.value : num [1:11, 1:11] 0.00 6.11e-10 9.38e-10 1.79e-07 1.78e-05 ...
##  $ lower.ci: num [1:11, 1:11] 1 -0.926 -0.923 -0.885 0.436 ...
##  $ upper.ci: num [1:11, 1:11] 1 -0.716 -0.708 -0.586 0.832 ...
##  - attr(*, "class")= chr "correlation"

使用其它相关系数算法。

## spearman
correlate(mtcars, cor.test = TRUE, method = "spearman") %>% str()
## List of 4
##  $ r       : num [1:11, 1:11] 1 -0.911 -0.909 -0.895 0.651 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ p.value : num [1:11, 1:11] 0.00 4.69e-13 6.37e-13 5.09e-12 5.38e-05 ...
##  $ lower.ci: NULL
##  $ upper.ci: NULL
##  - attr(*, "class")= chr "correlation"
## kendall
correlate(mtcars, cor.test = TRUE, method = "kendall") %>% str()
## List of 4
##  $ r       : num [1:11, 1:11] 1 -0.795 -0.768 -0.743 0.465 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ p.value : num [1:11, 1:11] 1.95e-15 2.25e-08 1.01e-09 4.33e-09 2.38e-04 ...
##  $ lower.ci: NULL
##  $ upper.ci: NULL
##  - attr(*, "class")= chr "correlation"

两个数据集的交叉相关系数计算。

library(vegan)
data("varespec")
data("varechem")
correlate(varechem, varespec) %>% str()
## List of 4
##  $ r       : num [1:14, 1:44] 0.19478 -0.08983 0.50892 0.00432 0.11448 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:14] "N" "P" "K" "Ca" ...
##   .. ..$ : chr [1:44] "Callvulg" "Empenigr" "Rhodtome" "Vaccmyrt" ...
##  $ p.value : NULL
##  $ lower.ci: NULL
##  $ upper.ci: NULL
##  - attr(*, "class")= chr "correlation"
correlate(varechem, varespec, cor.test = TRUE) %>% str()
## List of 4
##  $ r       : num [1:14, 1:44] 0.19478 -0.08983 0.50892 0.00432 0.11448 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:14] "N" "P" "K" "Ca" ...
##   .. ..$ : chr [1:44] "Callvulg" "Empenigr" "Rhodtome" "Vaccmyrt" ...
##  $ p.value : num [1:14, 1:44] 0.3617 0.6764 0.0111 0.984 0.5943 ...
##  $ lower.ci: num [1:14, 1:44] -0.226 -0.476 0.133 -0.4 -0.303 ...
##  $ upper.ci: num [1:14, 1:44] 0.555 0.325 0.757 0.407 0.495 ...
##  - attr(*, "class")= chr "correlation"
correlate(varechem, varespec, cor.test = TRUE, method = "spearman") %>% 
  str()
## List of 4
##  $ r       : num [1:14, 1:44] -0.2261 -0.0742 0.1921 -0.021 -0.0643 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:14] "N" "P" "K" "Ca" ...
##   .. ..$ : chr [1:44] "Callvulg" "Empenigr" "Rhodtome" "Vaccmyrt" ...
##  $ p.value : num [1:14, 1:44] 0.288 0.731 0.368 0.922 0.765 ...
##  $ lower.ci: NULL
##  $ upper.ci: NULL
##  - attr(*, "class")= chr "correlation"

有两个被问的最多的问题在这里说明:

  • ggcor的相关系数检验P值是否进行了校正?明确的说,没有,需要校正的请使用其它包。

  • correlate()计算效率如何?超过500个变量基本要蹦,我测试是500个变量,100个样本,需要40秒。

fast_correlate()函数

fast_correlate()函数是WGCNA::corAndPvalue()函数的简单封装,计算效率高,但是只支持pearson相关系数。ggcor不依赖WGCNA包,要使用该函数请先安装WGCNA包。

args(fast_correlate)
## function (x, y = NULL, use = "everything", ...) 
## NULL
fast_correlate(varechem, varespec) %>% str()
## List of 2
##  $ r      : num [1:14, 1:44] 0.19478 -0.08983 0.50892 0.00432 0.11448 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:14] "N" "P" "K" "Ca" ...
##   .. ..$ : chr [1:44] "Callvulg" "Empenigr" "Rhodtome" "Vaccmyrt" ...
##  $ p.value: num [1:14, 1:44] 0.3617 0.6764 0.0111 0.984 0.5943 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:14] "N" "P" "K" "Ca" ...
##   .. ..$ : chr [1:44] "Callvulg" "Empenigr" "Rhodtome" "Vaccmyrt" ...
##  - attr(*, "class")= chr "correlation"

常见包中的相关系数计算

有一些包对相关系数计算和检验做了专门的优化,更符合某些领域的专业需求,这里也简单的提及一下。

psych

psych包在心理学和医学领域的使用应该很广泛,提供了相关系数检验校正的一系列方法。理论上的内容我不清楚,大概就是知道有这些东西。

library(psych)
corr.test(mtcars) %>% str()
## List of 11
##  $ r     : num [1:11, 1:11] 1 -0.852 -0.848 -0.776 0.681 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ n     : num 32
##  $ t     : num [1:11, 1:11] Inf -8.92 -8.75 -6.74 5.1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ p     : num [1:11, 1:11] 0.00 6.11e-10 9.38e-10 1.79e-07 1.78e-05 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ se    : num [1:11, 1:11] 0 0.0955 0.0969 0.1151 0.1337 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ sef   : num 0.186
##  $ adjust: chr "holm"
##  $ sym   : logi TRUE
##  $ ci    :'data.frame':  55 obs. of  4 variables:
##   ..$ lower: num [1:55] -0.926 -0.923 -0.885 0.436 -0.934 ...
##   ..$ r    : num [1:55] -0.852 -0.848 -0.776 0.681 -0.868 ...
##   ..$ upper: num [1:55] -0.716 -0.708 -0.586 0.832 -0.744 ...
##   ..$ p    : num [1:55] 6.11e-10 9.38e-10 1.79e-07 1.78e-05 1.29e-10 ...
##  $ ci.adj:'data.frame':  55 obs. of  2 variables:
##   ..$ lower.adj: num [1:55] -0.954 -0.953 -0.928 0.238 -0.959 ...
##   ..$ upper.adj: num [1:55] -0.572 -0.562 -0.405 0.89 -0.61 ...
##  $ Call  : language corr.test(x = mtcars)
##  - attr(*, "class")= chr [1:2] "psych" "corr.test"

Hmisc

rcorr()函数要求输入必须为矩阵,所以要把数据框转化成matrix类。

library(Hmisc)
rcorr(data.matrix(mtcars)) %>% str()
## List of 3
##  $ r: num [1:11, 1:11] 1 -0.852 -0.848 -0.776 0.681 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ n: int [1:11, 1:11] 32 32 32 32 32 32 32 32 32 32 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ P: num [1:11, 1:11] NA 6.11e-10 9.38e-10 1.79e-07 1.78e-05 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  - attr(*, "class")= chr "rcorr"

WGCNA

这个包可以使用fast_correlate()函数,也可以直接使用corAndPvalue()函数。

library(WGCNA)
corAndPvalue(mtcars) %>% str()
## List of 5
##  $ cor : num [1:11, 1:11] 1 -0.852 -0.848 -0.776 0.681 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ p   : num [1:11, 1:11] 0.00 6.11e-10 9.38e-10 1.79e-07 1.78e-05 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ Z   : num [1:11, 1:11] Inf -6.92 -6.83 -5.67 4.55 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ t   : num [1:11, 1:11] Inf 8.92 8.75 6.74 5.1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##  $ nObs: num [1:11, 1:11] 32 32 32 32 32 32 32 32 32 32 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:11] "mpg" "cyl" "disp" "hp" ...

小结

这篇文章感觉和ggcor关联不大,但是是后续很多问题的基础,理解相关系数计算和检验的来龙去脉,能更好地选择合适的方法,也更容易理解,ggcor底层干了什么。

下期预告:数据转化。