厚缊

诹图——ggcor简介(二)

厚缊 / 2020-02-01


上一篇讲了ggcor默认是如何计算相关性的,也大概介绍了其它几个扩展包相关计算及检验函数,今天主要介绍如何把这些函数计算的结果转化成cor_tbl对象(即ggcor相关系数可视化的数据格式)。

基础函数

ggcor的可视化都是依赖ggplot2包,所以需要把相关系数矩阵转化成长数据(即data.frame对象),ggcor中完成这一工作的默认函数是cor_tbl()函数。前两个参数分别是相关系数矩阵和P值矩阵,extra.mat参数提供了一个接口,可以同时增加一些和相关系数矩阵维度一样的其它数据。

library(ggcor)
args(cor_tbl)
## function (corr, p.value = NULL, extra.mat = list(), type = "full", 
##     show.diag = TRUE, row.names = NULL, col.names = NULL, cluster = FALSE, 
##     ...) 
## NULL
cor(mtcars) %>% cor_tbl()
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.852      10       1
##  3 disp       mpg        -0.848       9       1
##  4 hp         mpg        -0.776       8       1
##  5 drat       mpg         0.681       7       1
##  6 wt         mpg        -0.868       6       1
##  7 qsec       mpg         0.419       5       1
##  8 vs         mpg         0.664       4       1
##  9 am         mpg         0.600       3       1
## 10 gear       mpg         0.480       2       1
## # … with 111 more rows

只保留上/下三角数据,其它行删除。type = "upper" 是保留上三角,默认情况下,保留了对角线,若要移除对角线上的数据,需设置show.diag = FALSE

cor(mtcars) %>% 
  cor_tbl(type = "upper") 
## # A tibble: 66 x 5
##    .row.names .col.names      r .row.id .col.id
##    <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 mpg        cyl        -0.852      11       2
##  3 cyl        cyl         1          10       2
##  4 mpg        disp       -0.848      11       3
##  5 cyl        disp        0.902      10       3
##  6 disp       disp        1           9       3
##  7 mpg        hp         -0.776      11       4
##  8 cyl        hp          0.832      10       4
##  9 disp       hp          0.791       9       4
## 10 hp         hp          1           8       4
## # … with 56 more rows
cor(mtcars) %>% 
  cor_tbl(type = "lower", show.diag = FALSE)
## # A tibble: 55 x 5
##    .row.names .col.names      r .row.id .col.id
##    <chr>      <chr>       <dbl>   <int>   <int>
##  1 cyl        mpg        -0.852      10       1
##  2 disp       mpg        -0.848       9       1
##  3 hp         mpg        -0.776       8       1
##  4 drat       mpg         0.681       7       1
##  5 wt         mpg        -0.868       6       1
##  6 qsec       mpg         0.419       5       1
##  7 vs         mpg         0.664       4       1
##  8 am         mpg         0.600       3       1
##  9 gear       mpg         0.480       2       1
## 10 carb       mpg        -0.551       1       1
## # … with 45 more rows

若需要对相关系数矩阵进行聚类,可设置cluster = TRUE最早公开的版本是支持分别按照行或者列聚类,而且同时支持对称或者不对称矩阵聚类,新版仅支持对对称的相关系数矩阵聚类,这样做的目的主要是为了和corrplot包保持一致性。

cor(mtcars) %>% 
  cor_tbl(cluster = TRUE)
## # A tibble: 121 x 5
##    .row.names .col.names       r .row.id .col.id
##  * <chr>      <chr>        <dbl>   <int>   <int>
##  1 carb       carb        1           11       1
##  2 wt         carb        0.428       10       1
##  3 hp         carb        0.750        9       1
##  4 cyl        carb        0.527        8       1
##  5 disp       carb        0.395        7       1
##  6 qsec       carb       -0.656        6       1
##  7 vs         carb       -0.570        5       1
##  8 mpg        carb       -0.551        4       1
##  9 drat       carb       -0.0908       3       1
## 10 am         carb        0.0575       2       1
## # … with 111 more rows

as_cor_tbl()泛型函数

as_cor_tbl()是一个S3的泛型函数,目前主要提供了六种方法。除了mantel_tbl对象外,其它几种均会调用cor_tbl()函数,也就是说cor_tbl()函数有的参数在这五种情况下依然支持。

methods(as_cor_tbl)
## [1] as_cor_tbl.corr.test*   as_cor_tbl.correlation* as_cor_tbl.data.frame* 
## [4] as_cor_tbl.default*     as_cor_tbl.mantel_tbl*  as_cor_tbl.matrix*     
## [7] as_cor_tbl.rcorr*      
## see '?methods' for accessing help and source code

correlation对象转化

correlation对象是ggcorcorrelate()fast_correlate()函数的输出结果,包含相关系数矩阵、P值矩阵、统计置信区间(针对前者),其结果能直接传递给as_cor_tbl()

correlate(mtcars) %>% as_cor_tbl() # 默认
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.852      10       1
##  3 disp       mpg        -0.848       9       1
##  4 hp         mpg        -0.776       8       1
##  5 drat       mpg         0.681       7       1
##  6 wt         mpg        -0.868       6       1
##  7 qsec       mpg         0.419       5       1
##  8 vs         mpg         0.664       4       1
##  9 am         mpg         0.600       3       1
## 10 gear       mpg         0.480       2       1
## # … with 111 more rows
correlate(mtcars) %>% as_cor_tbl(type = "upper") # 上三角
## # A tibble: 66 x 5
##    .row.names .col.names      r .row.id .col.id
##    <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 mpg        cyl        -0.852      11       2
##  3 cyl        cyl         1          10       2
##  4 mpg        disp       -0.848      11       3
##  5 cyl        disp        0.902      10       3
##  6 disp       disp        1           9       3
##  7 mpg        hp         -0.776      11       4
##  8 cyl        hp          0.832      10       4
##  9 disp       hp          0.791       9       4
## 10 hp         hp          1           8       4
## # … with 56 more rows
correlate(mtcars, cor.test = TRUE) %>% 
  as_cor_tbl(type = "lower", show.diag = FALSE) # 下三角不带对角线
## # A tibble: 55 x 8
##    .row.names .col.names      r  p.value upper.ci lower.ci .row.id .col.id
##    <chr>      <chr>       <dbl>    <dbl>    <dbl>    <dbl>   <int>   <int>
##  1 cyl        mpg        -0.852 6.11e-10   -0.716  -0.926       10       1
##  2 disp       mpg        -0.848 9.38e-10   -0.708  -0.923        9       1
##  3 hp         mpg        -0.776 1.79e- 7   -0.586  -0.885        8       1
##  4 drat       mpg         0.681 1.78e- 5    0.832   0.436        7       1
##  5 wt         mpg        -0.868 1.29e-10   -0.744  -0.934        6       1
##  6 qsec       mpg         0.419 1.71e- 2    0.670   0.0820       5       1
##  7 vs         mpg         0.664 3.42e- 5    0.822   0.410        4       1
##  8 am         mpg         0.600 2.85e- 4    0.784   0.318        3       1
##  9 gear       mpg         0.480 5.40e- 3    0.710   0.158        2       1
## 10 carb       mpg        -0.551 1.08e- 3   -0.250  -0.755        1       1
## # … with 45 more rows
library(vegan)
data("varespec")
data("varechem")
fast_correlate(varechem, varespec) %>% 
  as_cor_tbl() # 非对称相关系数矩阵
## # A tibble: 616 x 6
##    .row.names .col.names        r p.value .row.id .col.id
##  * <chr>      <chr>         <dbl>   <dbl>   <int>   <int>
##  1 N          Callvulg    0.195    0.362       14       1
##  2 P          Callvulg   -0.0898   0.676       13       1
##  3 K          Callvulg    0.509    0.0111      12       1
##  4 Ca         Callvulg    0.00432  0.984       11       1
##  5 Mg         Callvulg    0.114    0.594       10       1
##  6 S          Callvulg    0.401    0.0520       9       1
##  7 Al         Callvulg    0.330    0.116        8       1
##  8 Fe         Callvulg    0.168    0.433        7       1
##  9 Mn         Callvulg    0.0354   0.870        6       1
## 10 Zn         Callvulg    0.0902   0.675        5       1
## # … with 606 more rows

corr.test对象转化

有时候需要校正相关系数检验,psych包中的corr.test()函数是不错的选择,其返回结果就是corr.test对象,能直接调用as_cor_tbl()函数进行转化。

library(psych) # 需要就按照这个包,不需要可以跳过
corr.test(mtcars) %>% as_cor_tbl() # 默认
## # A tibble: 121 x 6
##    .row.names .col.names      r  p.value .row.id .col.id
##  * <chr>      <chr>       <dbl>    <dbl>   <int>   <int>
##  1 mpg        mpg         1     0.            11       1
##  2 cyl        mpg        -0.852 6.11e-10      10       1
##  3 disp       mpg        -0.848 9.38e-10       9       1
##  4 hp         mpg        -0.776 1.79e- 7       8       1
##  5 drat       mpg         0.681 1.78e- 5       7       1
##  6 wt         mpg        -0.868 1.29e-10       6       1
##  7 qsec       mpg         0.419 1.71e- 2       5       1
##  8 vs         mpg         0.664 3.42e- 5       4       1
##  9 am         mpg         0.600 2.85e- 4       3       1
## 10 gear       mpg         0.480 5.40e- 3       2       1
## # … with 111 more rows
corr.test(mtcars) %>% as_cor_tbl(type = "upper", cluster = TRUE)
## # A tibble: 66 x 6
##    .row.names .col.names     r       p.value .row.id .col.id
##    <chr>      <chr>      <dbl>         <dbl>   <int>   <int>
##  1 carb       carb       1     0                  11       1
##  2 carb       wt         0.428 0.0146             11       2
##  3 wt         wt         1     0                  10       2
##  4 carb       hp         0.750 0.000000783        11       3
##  5 wt         hp         0.659 0.0000415          10       3
##  6 hp         hp         1     0                   9       3
##  7 carb       cyl        0.527 0.00194            11       4
##  8 wt         cyl        0.782 0.000000122        10       4
##  9 hp         cyl        0.832 0.00000000348       9       4
## 10 cyl        cyl        1     0                   8       4
## # … with 56 more rows

rcorr对象转化

Hmisc中的rcorr()函数也能计算相关系数,而且比R自带的函数效率更高,其结果是rcorr对象,也是as_cor_tbl()直接支持的类型。

library(Hmisc) # 需要就按照这个包,不需要可以跳过
rcorr(data.matrix(mtcars)) %>% as_cor_tbl() # 默认
## # A tibble: 121 x 6
##    .row.names .col.names      r  p.value .row.id .col.id
##  * <chr>      <chr>       <dbl>    <dbl>   <int>   <int>
##  1 mpg        mpg         1     0.            11       1
##  2 cyl        mpg        -0.852 6.11e-10      10       1
##  3 disp       mpg        -0.848 9.38e-10       9       1
##  4 hp         mpg        -0.776 1.79e- 7       8       1
##  5 drat       mpg         0.681 1.78e- 5       7       1
##  6 wt         mpg        -0.868 1.29e-10       6       1
##  7 qsec       mpg         0.419 1.71e- 2       5       1
##  8 vs         mpg         0.664 3.42e- 5       4       1
##  9 am         mpg         0.600 2.85e- 4       3       1
## 10 gear       mpg         0.480 5.40e- 3       2       1
## # … with 111 more rows
rcorr(data.matrix(mtcars)) %>% 
  as_cor_tbl(type = "lower", show.diag = FALSE)
## # A tibble: 55 x 6
##    .row.names .col.names      r  p.value .row.id .col.id
##    <chr>      <chr>       <dbl>    <dbl>   <int>   <int>
##  1 cyl        mpg        -0.852 6.11e-10      10       1
##  2 disp       mpg        -0.848 9.38e-10       9       1
##  3 hp         mpg        -0.776 1.79e- 7       8       1
##  4 drat       mpg         0.681 1.78e- 5       7       1
##  5 wt         mpg        -0.868 1.29e-10       6       1
##  6 qsec       mpg         0.419 1.71e- 2       5       1
##  7 vs         mpg         0.664 3.42e- 5       4       1
##  8 am         mpg         0.600 2.85e- 4       3       1
##  9 gear       mpg         0.480 5.40e- 3       2       1
## 10 carb       mpg        -0.551 1.08e- 3       1       1
## # … with 45 more rows

测试下correlate()函数(基于R自带函数)和rcorr()的计算效率可以发现,rcorr()差不多快了400倍。

library(microbenchmark)
m <- matrix(rnorm(100*100), nrow = 100)
microbenchmark(
  corr1 = correlate(m, cor.test = TRUE),
  corr2 = rcorr(m),
  times = 10
)
## Unit: milliseconds
##   expr         min          lq        mean      median          uq
##  corr1 1531.417954 1554.849422 1616.620893 1567.063541 1659.959056
##  corr2    4.121187    4.546143    4.567405    4.617854    4.670538
##          max neval
##  1812.434420    10
##     4.679552    10

matrixdata.frame对象转化

对于matrixdata.frame对象,as_cor_tbl()直接当成相关系数矩阵处理,若是原始数据,切不可直接调用

mat <- cor(mtcars)
df <- as.data.frame(mat)
as_cor_tbl(mat)
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.852      10       1
##  3 disp       mpg        -0.848       9       1
##  4 hp         mpg        -0.776       8       1
##  5 drat       mpg         0.681       7       1
##  6 wt         mpg        -0.868       6       1
##  7 qsec       mpg         0.419       5       1
##  8 vs         mpg         0.664       4       1
##  9 am         mpg         0.600       3       1
## 10 gear       mpg         0.480       2       1
## # … with 111 more rows
as_cor_tbl(df)
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.852      10       1
##  3 disp       mpg        -0.848       9       1
##  4 hp         mpg        -0.776       8       1
##  5 drat       mpg         0.681       7       1
##  6 wt         mpg        -0.868       6       1
##  7 qsec       mpg         0.419       5       1
##  8 vs         mpg         0.664       4       1
##  9 am         mpg         0.600       3       1
## 10 gear       mpg         0.480       2       1
## # … with 111 more rows

fortify_cor()函数

一般情况

fortify_cor()as_cor_tbl()函数的深度包装,能根据输入对象类型灵活的进行数据转换。对于as_cor_tbl()支持的类型,直接调用as_cor_tbl()函数,对于原始数据,先调用correlate()函数进行相关系数计算和检验,然后调用as_cor_tbl()函数转换。

fortify_cor(mtcars)
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.852      10       1
##  3 disp       mpg        -0.848       9       1
##  4 hp         mpg        -0.776       8       1
##  5 drat       mpg         0.681       7       1
##  6 wt         mpg        -0.868       6       1
##  7 qsec       mpg         0.419       5       1
##  8 vs         mpg         0.664       4       1
##  9 am         mpg         0.600       3       1
## 10 gear       mpg         0.480       2       1
## # … with 111 more rows
correlate(mtcars) %>% fortify_cor()
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.852      10       1
##  3 disp       mpg        -0.848       9       1
##  4 hp         mpg        -0.776       8       1
##  5 drat       mpg         0.681       7       1
##  6 wt         mpg        -0.868       6       1
##  7 qsec       mpg         0.419       5       1
##  8 vs         mpg         0.664       4       1
##  9 am         mpg         0.600       3       1
## 10 gear       mpg         0.480       2       1
## # … with 111 more rows
rcorr(data.matrix(mtcars)) %>% fortify_cor()
## # A tibble: 121 x 6
##    .row.names .col.names      r  p.value .row.id .col.id
##  * <chr>      <chr>       <dbl>    <dbl>   <int>   <int>
##  1 mpg        mpg         1     0.            11       1
##  2 cyl        mpg        -0.852 6.11e-10      10       1
##  3 disp       mpg        -0.848 9.38e-10       9       1
##  4 hp         mpg        -0.776 1.79e- 7       8       1
##  5 drat       mpg         0.681 1.78e- 5       7       1
##  6 wt         mpg        -0.868 1.29e-10       6       1
##  7 qsec       mpg         0.419 1.71e- 2       5       1
##  8 vs         mpg         0.664 3.42e- 5       4       1
##  9 am         mpg         0.600 2.85e- 4       3       1
## 10 gear       mpg         0.480 5.40e- 3       2       1
## # … with 111 more rows
corr.test(mtcars) %>% fortify_cor()
## # A tibble: 121 x 6
##    .row.names .col.names      r  p.value .row.id .col.id
##  * <chr>      <chr>       <dbl>    <dbl>   <int>   <int>
##  1 mpg        mpg         1     0.            11       1
##  2 cyl        mpg        -0.852 6.11e-10      10       1
##  3 disp       mpg        -0.848 9.38e-10       9       1
##  4 hp         mpg        -0.776 1.79e- 7       8       1
##  5 drat       mpg         0.681 1.78e- 5       7       1
##  6 wt         mpg        -0.868 1.29e-10       6       1
##  7 qsec       mpg         0.419 1.71e- 2       5       1
##  8 vs         mpg         0.664 3.42e- 5       4       1
##  9 am         mpg         0.600 2.85e- 4       3       1
## 10 gear       mpg         0.480 5.40e- 3       2       1
## # … with 111 more rows

若传入的是相关系数矩阵,需要设置is.cor = TRUE,否则会误当原始数据处理。

corr <- cor(mtcars)
fortify_cor(corr, is.cor = TRUE) # 正确
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.852      10       1
##  3 disp       mpg        -0.848       9       1
##  4 hp         mpg        -0.776       8       1
##  5 drat       mpg         0.681       7       1
##  6 wt         mpg        -0.868       6       1
##  7 qsec       mpg         0.419       5       1
##  8 vs         mpg         0.664       4       1
##  9 am         mpg         0.600       3       1
## 10 gear       mpg         0.480       2       1
## # … with 111 more rows
fortify_cor(corr) # 错误
## # A tibble: 121 x 5
##    .row.names .col.names      r .row.id .col.id
##  * <chr>      <chr>       <dbl>   <int>   <int>
##  1 mpg        mpg         1          11       1
##  2 cyl        mpg        -0.991      10       1
##  3 disp       mpg        -0.993       9       1
##  4 hp         mpg        -0.956       8       1
##  5 drat       mpg         0.939       7       1
##  6 wt         mpg        -0.986       6       1
##  7 qsec       mpg         0.708       5       1
##  8 vs         mpg         0.932       4       1
##  9 am         mpg         0.827       3       1
## 10 gear       mpg         0.768       2       1
## # … with 111 more rows

样本分组

这部分功能是测试性的,主要考虑的问题是如何处理需要针对不同的样本进行相关系数计算和检验的情况。fortify_cor()有一个group参数,可以传入一个样本(即原始数据行)分组索引,然后根据索引信息分别计算和转化。

group <- iris[[5]]
fortify_cor(iris[-5], group = group)
## # A tibble: 48 x 6
##    .row.names   .col.names       r .row.id .col.id .group
##    <chr>        <chr>        <dbl>   <int>   <int> <chr> 
##  1 Sepal.Length Sepal.Length 1           4       1 setosa
##  2 Sepal.Width  Sepal.Length 0.743       3       1 setosa
##  3 Petal.Length Sepal.Length 0.267       2       1 setosa
##  4 Petal.Width  Sepal.Length 0.278       1       1 setosa
##  5 Sepal.Length Sepal.Width  0.743       4       2 setosa
##  6 Sepal.Width  Sepal.Width  1           3       2 setosa
##  7 Petal.Length Sepal.Width  0.178       2       2 setosa
##  8 Petal.Width  Sepal.Width  0.233       1       2 setosa
##  9 Sepal.Length Petal.Length 0.267       4       3 setosa
## 10 Sepal.Width  Petal.Length 0.178       3       3 setosa
## # … with 38 more rows
fortify_cor(iris[-5], group = group, type = "upper")
## # A tibble: 18 x 6
##    .row.names   .col.names       r .row.id .col.id .group    
##    <chr>        <chr>        <dbl>   <int>   <int> <chr>     
##  1 Sepal.Length Sepal.Width  0.743       4       2 setosa    
##  2 Sepal.Length Petal.Length 0.267       4       3 setosa    
##  3 Sepal.Width  Petal.Length 0.178       3       3 setosa    
##  4 Sepal.Length Petal.Width  0.278       4       4 setosa    
##  5 Sepal.Width  Petal.Width  0.233       3       4 setosa    
##  6 Petal.Length Petal.Width  0.332       2       4 setosa    
##  7 Sepal.Length Sepal.Width  0.526       4       2 versicolor
##  8 Sepal.Length Petal.Length 0.754       4       3 versicolor
##  9 Sepal.Width  Petal.Length 0.561       3       3 versicolor
## 10 Sepal.Length Petal.Width  0.546       4       4 versicolor
## 11 Sepal.Width  Petal.Width  0.664       3       4 versicolor
## 12 Petal.Length Petal.Width  0.787       2       4 versicolor
## 13 Sepal.Length Sepal.Width  0.457       4       2 virginica 
## 14 Sepal.Length Petal.Length 0.864       4       3 virginica 
## 15 Sepal.Width  Petal.Length 0.401       3       3 virginica 
## 16 Sepal.Length Petal.Width  0.281       4       4 virginica 
## 17 Sepal.Width  Petal.Width  0.538       3       4 virginica 
## 18 Petal.Length Petal.Width  0.322       2       4 virginica

小结

ggcor相关系数可视化完全依赖于cor_tbl对象,也就是依赖上述函数把其它对象转成cor_tbl,因此理解上述函数是灵活使用ggcor的基础。

下期预告:相关系数矩阵热图。