厚缊

诹图——ggcor新版本变化

厚缊 / 2019-12-22


ggcor从酝酿到现在,真是骑虎难下,一言难尽。由于前期是临时起意做这个包,所有的代码都没有经过认真的思考和严格的逻辑设计,堆砌的情况非常严重。在嘉杰的帮助下,重新梳理了下逻辑结构,整体重写了几个重要函数的代码,现在无论从代码的可读性和效率来说,应该都有不少提升。

devtools::install_github("houyunhuang/ggcor")

变化(一)

打包了相关性计算和检验的correlate()函数,ggcor中对原始数据均通过此函数来进行计算和检验。cor.test参数是是否进行相关性检验的逻辑值,其它参数均和cor()cor.test()中一致。会根据条件返回r、p.value、lower.ci、upper.ci矩阵。

library(ggcor)
args(correlate)
## function (x, y = NULL, cor.test = FALSE, method = "pearson", 
##     use = "everything", ...) 
## NULL
cor <- correlate(mtcars, cor.test = TRUE)
names(cor)
## [1] "r"        "p.value"  "lower.ci" "upper.ci"

变化(二)

前面的几篇测试性的文章就已经反复说过,ggcor保持巨大灵活性的基础是as_cor_tbl()函数,能打通和其它包之间的联系。

matrixdata.frame类的转换

我们都知道,cor()对相关性进行计算后返回相关系数矩阵,那么最基础的形式也就是如何把相关系数矩阵转换为cor_tbl对象。

1.通过cor()计算mtcars数据集的相关性矩阵;

2.用as_cor_tbl()转换成cor_tbl对象。

m <- cor(mtcars)
as_cor_tbl(m)
## # A tibble: 121 x 5
##    idx   idy        r     x     y
##  * <chr> <chr>  <dbl> <int> <int>
##  1 mpg   mpg    1         1    11
##  2 mpg   cyl   -0.852     1    10
##  3 mpg   disp  -0.848     1     9
##  4 mpg   hp    -0.776     1     8
##  5 mpg   drat   0.681     1     7
##  6 mpg   wt    -0.868     1     6
##  7 mpg   qsec   0.419     1     5
##  8 mpg   vs     0.664     1     4
##  9 mpg   am     0.600     1     3
## 10 mpg   gear   0.480     1     2
## # … with 111 more rows

若是用过ggcor0.8.0之前的版本的,可能会发现新版本数据多了两列idx、idy两个变量,直接列出了每一行是哪两个变量之间的相关性,对于用户来说更友好。当然,这样的变化还有个考虑是方便正在增加的网络分析。

这个时候可能有疑问,若是还要增加p value变量怎么整?方法很简单,看下面的例子。当然,还可以类似的增加lower.ci、upper.ci变量。

r <- cor$r
p.value <- cor$p.value
as_cor_tbl(r, p.value = p.value)
## # A tibble: 121 x 6
##    idx   idy        r  p.value     x     y
##  * <chr> <chr>  <dbl>    <dbl> <int> <int>
##  1 mpg   mpg    1     0.           1    11
##  2 mpg   cyl   -0.852 6.11e-10     1    10
##  3 mpg   disp  -0.848 9.38e-10     1     9
##  4 mpg   hp    -0.776 1.79e- 7     1     8
##  5 mpg   drat   0.681 1.78e- 5     1     7
##  6 mpg   wt    -0.868 1.29e-10     1     6
##  7 mpg   qsec   0.419 1.71e- 2     1     5
##  8 mpg   vs     0.664 3.42e- 5     1     4
##  9 mpg   am     0.600 2.85e- 4     1     3
## 10 mpg   gear   0.480 5.40e- 3     1     2
## # … with 111 more rows

尽管我们从R中直接获得相关系数矩阵为data.frame对象的可能性很小,但是有一种可能是你用其它软件计算了相关系数矩阵,然后用read.csv()之类的函数读进R来做可视化,这时候就比较有用。

r_df <- as.data.frame(cor$r) ## 把矩阵转成df
as_cor_tbl(r_df)
## # A tibble: 121 x 5
##    idx   idy        r     x     y
##  * <chr> <chr>  <dbl> <int> <int>
##  1 mpg   mpg    1         1    11
##  2 mpg   cyl   -0.852     1    10
##  3 mpg   disp  -0.848     1     9
##  4 mpg   hp    -0.776     1     8
##  5 mpg   drat   0.681     1     7
##  6 mpg   wt    -0.868     1     6
##  7 mpg   qsec   0.419     1     5
##  8 mpg   vs     0.664     1     4
##  9 mpg   am     0.600     1     3
## 10 mpg   gear   0.480     1     2
## # … with 111 more rows

其它类数据转换

看了前面的转换,肯定会觉得as_cor_tbl()有的呆滞,还要分别指定每个变量矩阵。其实从本质上说,上面的matrixdata.frame方法仅仅是基础,一般情况下我们都比较少直接使用。

correlate对象转换

as_cor_tbl(cor)
## # A tibble: 121 x 8
##    idx   idy        r  p.value lower.ci upper.ci     x     y
##  * <chr> <chr>  <dbl>    <dbl>    <dbl>    <dbl> <int> <int>
##  1 mpg   mpg    1     0.         1         1         1    11
##  2 mpg   cyl   -0.852 6.11e-10  -0.926    -0.716     1    10
##  3 mpg   disp  -0.848 9.38e-10  -0.923    -0.708     1     9
##  4 mpg   hp    -0.776 1.79e- 7  -0.885    -0.586     1     8
##  5 mpg   drat   0.681 1.78e- 5   0.436     0.832     1     7
##  6 mpg   wt    -0.868 1.29e-10  -0.934    -0.744     1     6
##  7 mpg   qsec   0.419 1.71e- 2   0.0820    0.670     1     5
##  8 mpg   vs     0.664 3.42e- 5   0.410     0.822     1     4
##  9 mpg   am     0.600 2.85e- 4   0.318     0.784     1     3
## 10 mpg   gear   0.480 5.40e- 3   0.158     0.710     1     2
## # … with 111 more rows

rcorr对象转换

rcorr对象是Hmisc包中rcorr()函数返回的结果,其中包含了相关系数矩阵r和统计检验p value。

library(Hmisc)
rcorr <- rcorr(data.matrix(mtcars)) ## 有点变态,竟然必须是矩阵输入
as_cor_tbl(rcorr)
## # A tibble: 121 x 6
##    idx   idy        r  p.value     x     y
##  * <chr> <chr>  <dbl>    <dbl> <int> <int>
##  1 mpg   mpg    1     0.           1    11
##  2 mpg   cyl   -0.852 6.11e-10     1    10
##  3 mpg   disp  -0.848 9.38e-10     1     9
##  4 mpg   hp    -0.776 1.79e- 7     1     8
##  5 mpg   drat   0.681 1.78e- 5     1     7
##  6 mpg   wt    -0.868 1.29e-10     1     6
##  7 mpg   qsec   0.419 1.71e- 2     1     5
##  8 mpg   vs     0.664 3.42e- 5     1     4
##  9 mpg   am     0.600 2.85e- 4     1     3
## 10 mpg   gear   0.480 5.40e- 3     1     2
## # … with 111 more rows

corr.test对象转换

corr.test对象是psych包中corr.test函数返回的结果。

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

list对象的转换

可能用的最多的是WGCNA包中corAndPvalue()函数,该函数对于相关性检验做了优化,计算非常快。

library(WGCNA)
cor_p <- corAndPvalue(mtcars)
as_cor_tbl(cor_p)
## # A tibble: 121 x 6
##    idx   idy        r  p.value     x     y
##  * <chr> <chr>  <dbl>    <dbl> <int> <int>
##  1 mpg   mpg    1     0.           1    11
##  2 mpg   cyl   -0.852 6.11e-10     1    10
##  3 mpg   disp  -0.848 9.38e-10     1     9
##  4 mpg   hp    -0.776 1.79e- 7     1     8
##  5 mpg   drat   0.681 1.78e- 5     1     7
##  6 mpg   wt    -0.868 1.29e-10     1     6
##  7 mpg   qsec   0.419 1.71e- 2     1     5
##  8 mpg   vs     0.664 3.42e- 5     1     4
##  9 mpg   am     0.600 2.85e- 4     1     3
## 10 mpg   gear   0.480 5.40e- 3     1     2
## # … with 111 more rows

变化(三)

为了保证延续性,新版本保留了fortify_cor()函数,并在里面做了一些功能上的扩展。

as_cor_tbl()函数的深度封装

fortify_cor(cor)
## # A tibble: 121 x 8
##    idx   idy        r  p.value lower.ci upper.ci     x     y
##  * <chr> <chr>  <dbl>    <dbl>    <dbl>    <dbl> <int> <int>
##  1 mpg   mpg    1     0.         1         1         1    11
##  2 mpg   cyl   -0.852 6.11e-10  -0.926    -0.716     1    10
##  3 mpg   disp  -0.848 9.38e-10  -0.923    -0.708     1     9
##  4 mpg   hp    -0.776 1.79e- 7  -0.885    -0.586     1     8
##  5 mpg   drat   0.681 1.78e- 5   0.436     0.832     1     7
##  6 mpg   wt    -0.868 1.29e-10  -0.934    -0.744     1     6
##  7 mpg   qsec   0.419 1.71e- 2   0.0820    0.670     1     5
##  8 mpg   vs     0.664 3.42e- 5   0.410     0.822     1     4
##  9 mpg   am     0.600 2.85e- 4   0.318     0.784     1     3
## 10 mpg   gear   0.480 5.40e- 3   0.158     0.710     1     2
## # … with 111 more rows
fortify_cor(rcorr)
## # A tibble: 121 x 6
##    idx   idy        r  p.value     x     y
##  * <chr> <chr>  <dbl>    <dbl> <int> <int>
##  1 mpg   mpg    1     0.           1    11
##  2 mpg   cyl   -0.852 6.11e-10     1    10
##  3 mpg   disp  -0.848 9.38e-10     1     9
##  4 mpg   hp    -0.776 1.79e- 7     1     8
##  5 mpg   drat   0.681 1.78e- 5     1     7
##  6 mpg   wt    -0.868 1.29e-10     1     6
##  7 mpg   qsec   0.419 1.71e- 2     1     5
##  8 mpg   vs     0.664 3.42e- 5     1     4
##  9 mpg   am     0.600 2.85e- 4     1     3
## 10 mpg   gear   0.480 5.40e- 3     1     2
## # … with 111 more rows
fortify_cor(corr.test)
## # A tibble: 121 x 6
##    idx   idy        r  p.value     x     y
##  * <chr> <chr>  <dbl>    <dbl> <int> <int>
##  1 mpg   mpg    1     0.           1    11
##  2 mpg   cyl   -0.852 6.11e-10     1    10
##  3 mpg   disp  -0.848 9.38e-10     1     9
##  4 mpg   hp    -0.776 1.79e- 7     1     8
##  5 mpg   drat   0.681 1.78e- 5     1     7
##  6 mpg   wt    -0.868 1.29e-10     1     6
##  7 mpg   qsec   0.419 1.71e- 2     1     5
##  8 mpg   vs     0.664 3.42e- 5     1     4
##  9 mpg   am     0.600 2.85e- 4     1     3
## 10 mpg   gear   0.480 5.40e- 3     1     2
## # … with 111 more rows
fortify_cor(cor_p)
## # A tibble: 3,025 x 5
##    idx     idy           r     x     y
##  * <chr>   <chr>     <dbl> <int> <int>
##  1 cor.mpg cor.mpg   1.000     1    55
##  2 cor.mpg cor.cyl  -0.991     1    54
##  3 cor.mpg cor.disp -0.993     1    53
##  4 cor.mpg cor.hp   -0.956     1    52
##  5 cor.mpg cor.drat  0.939     1    51
##  6 cor.mpg cor.wt   -0.986     1    50
##  7 cor.mpg cor.qsec  0.708     1    49
##  8 cor.mpg cor.vs    0.932     1    48
##  9 cor.mpg cor.am    0.827     1    47
## 10 cor.mpg cor.gear  0.768     1    46
## # … with 3,015 more rows

内部计算

这里需要注意,当fortify_cor()输入数据是matrix或者data.frame对象时,若是相关系数矩阵,,要设置is.cor = TRUE参数,否则当成原始数据处理。

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

单个数据集

fortify_cor(mtcars)
## # A tibble: 121 x 5
##    idx   idy        r     x     y
##  * <chr> <chr>  <dbl> <int> <int>
##  1 mpg   mpg    1         1    11
##  2 mpg   cyl   -0.852     1    10
##  3 mpg   disp  -0.848     1     9
##  4 mpg   hp    -0.776     1     8
##  5 mpg   drat   0.681     1     7
##  6 mpg   wt    -0.868     1     6
##  7 mpg   qsec   0.419     1     5
##  8 mpg   vs     0.664     1     4
##  9 mpg   am     0.600     1     3
## 10 mpg   gear   0.480     1     2
## # … with 111 more rows

两个数据集

library(vegan)
data("varechem")
data("varespec")
fortify_cor(varechem, varespec)
## # A tibble: 616 x 5
##    idx      idy          r     x     y
##  * <chr>    <chr>    <dbl> <int> <int>
##  1 Callvulg N      0.195       1    14
##  2 Callvulg P     -0.0898      1    13
##  3 Callvulg K      0.509       1    12
##  4 Callvulg Ca     0.00432     1    11
##  5 Callvulg Mg     0.114       1    10
##  6 Callvulg S      0.401       1     9
##  7 Callvulg Al     0.330       1     8
##  8 Callvulg Fe     0.168       1     7
##  9 Callvulg Mn     0.0354      1     6
## 10 Callvulg Zn     0.0902      1     5
## # … with 606 more rows

样本分组

我不清楚这个功能有多少用处,但是我一直觉得,对于不同的样本数据混在一起做相关性是没有任何意义的,所有就增加了这个鸡肋的功能。样本分组也就是将原始数据按照行分成子集,然后分别针对每个子集中的变量计算相关性。

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
fortify_cor(iris[1:4], group = iris$Species)
## # A tibble: 48 x 6
##    idx          idy              r     x     y group 
##    <chr>        <chr>        <dbl> <int> <int> <chr> 
##  1 Sepal.Length Sepal.Length 1         1     4 setosa
##  2 Sepal.Length Sepal.Width  0.743     1     3 setosa
##  3 Sepal.Length Petal.Length 0.267     1     2 setosa
##  4 Sepal.Length Petal.Width  0.278     1     1 setosa
##  5 Sepal.Width  Sepal.Length 0.743     2     4 setosa
##  6 Sepal.Width  Sepal.Width  1         2     3 setosa
##  7 Sepal.Width  Petal.Length 0.178     2     2 setosa
##  8 Sepal.Width  Petal.Width  0.233     2     1 setosa
##  9 Petal.Length Sepal.Length 0.267     3     4 setosa
## 10 Petal.Length Sepal.Width  0.178     3     3 setosa
## # … with 38 more rows

对于两个数据集的情况也是可以的。

fortify_cor(varechem, varespec, group = rep(LETTERS[1:2], 12))
## # A tibble: 1,232 x 6
##    idx      idy          r     x     y group
##    <chr>    <chr>    <dbl> <int> <int> <chr>
##  1 Callvulg N     -0.137       1    14 A    
##  2 Callvulg P     -0.203       1    13 A    
##  3 Callvulg K      0.132       1    12 A    
##  4 Callvulg Ca     0.265       1    11 A    
##  5 Callvulg Mg     0.0463      1    10 A    
##  6 Callvulg S     -0.00143     1     9 A    
##  7 Callvulg Al    -0.102       1     8 A    
##  8 Callvulg Fe    -0.214       1     7 A    
##  9 Callvulg Mn    -0.0151      1     6 A    
## 10 Callvulg Zn    -0.0156      1     5 A    
## # … with 1,222 more rows

变化(四)

还有一大块的内容就是聚类方法变化,原版本依赖seriation包增加了一堆没用的聚类功能,现在完全简化了。目前只能对对称性相关系数矩阵聚类,聚类方法是调用stats::hclust()函数,当然为了简化,封装成了matrix_order()函数。

as_cor_tbl(m, cluster = TRUE) # 聚类
## # A tibble: 121 x 5
##    idx   idy         r     x     y
##  * <chr> <chr>   <dbl> <int> <int>
##  1 carb  carb   1          1    11
##  2 carb  wt     0.428      1    10
##  3 carb  hp     0.750      1     9
##  4 carb  cyl    0.527      1     8
##  5 carb  disp   0.395      1     7
##  6 carb  qsec  -0.656      1     6
##  7 carb  vs    -0.570      1     5
##  8 carb  mpg   -0.551      1     4
##  9 carb  drat  -0.0908     1     3
## 10 carb  am     0.0575     1     2
## # … with 111 more rows
as_cor_tbl(m, cluster = TRUE, cluster.method = "ward.D2")
## # A tibble: 121 x 5
##    idx   idy        r     x     y
##  * <chr> <chr>  <dbl> <int> <int>
##  1 wt    wt     1         1    11
##  2 wt    cyl    0.782     1    10
##  3 wt    disp   0.888     1     9
##  4 wt    hp     0.659     1     8
##  5 wt    carb   0.428     1     7
##  6 wt    qsec  -0.175     1     6
##  7 wt    vs    -0.555     1     5
##  8 wt    mpg   -0.868     1     4
##  9 wt    drat  -0.712     1     3
## 10 wt    am    -0.692     1     2
## # … with 111 more rows
fortify_cor(varechem, varespec, cluster = TRUE) # 非对称,不支持聚类
## Warning: 'cluster' just spports for symmetric correlation matrix.
## # A tibble: 616 x 5
##    idx      idy          r     x     y
##  * <chr>    <chr>    <dbl> <int> <int>
##  1 Callvulg N      0.195       1    14
##  2 Callvulg P     -0.0898      1    13
##  3 Callvulg K      0.509       1    12
##  4 Callvulg Ca     0.00432     1    11
##  5 Callvulg Mg     0.114       1    10
##  6 Callvulg S      0.401       1     9
##  7 Callvulg Al     0.330       1     8
##  8 Callvulg Fe     0.168       1     7
##  9 Callvulg Mn     0.0354      1     6
## 10 Callvulg Zn     0.0902      1     5
## # … with 606 more rows

对于前面提到几个类型的数据转换和计算,均可以设置聚类参数。

变化(五)

ggcor所有的工作都是为可视化服务的,绘图函数的变化可能才会真正的影响使用。最大的变化原ggcor()函数变成了quickcor()函数,现在的ggcor()函数不会在内部处理数据(除了坐标轴相关),也不会处理颜色映射。这么做主要是基于两方面的考虑:一是在画图主函数里面塞进去太多的东西,导致参数膨胀,需要一个更清爽的函数让熟悉ggplot2的人用起来更自由,这就是现在的ggcor()函数,尽量保持了ggplot2原汁原味的风格;二是对于一般探索性分析,或者不熟悉R的人来说,更希望用最少的代码来出图,这也就是quickcor()函数的作用。(PS:之所有没有叫qggcor()或者quickplot()之类的函数,一个是命名冲突,还有一个是念起来拗口。)

quickcor()函数

主要参数保留了原版本ggcor()的参数设置方式,后面通过一些例子来讲解。

args(quickcor)
## function (x, y = NULL, mapping = NULL, fill.colours = NULL, fill.bin = FALSE, 
##     grid.colour = "grey50", grid.size = 0.25, axis.x.position = "auto", 
##     axis.y.position = "auto", axis.label.drop = TRUE, legend.title = "corr", 
##     legend.position = "auto", legend.breaks = NULL, legend.labels = NULL, 
##     ...) 
## NULL

初始化绘图函数

quickcor(mtcars)

quickcor(mtcars, type = "upper")

quickcor(mtcars, type = "lower", show.diag = FALSE)

主要图层函数

quickcor(mtcars) + geom_colour()

quickcor(mtcars, type = "upper") + geom_circle2()

quickcor(mtcars, type = "lower", show.diag = FALSE) + geom_ellipse2()

quickcor(mtcars, cluster = TRUE) + geom_square()

quickcor(mtcars, cor.test = TRUE) + geom_confbox()

quickcor(mtcars, cor.test = TRUE) + geom_colour() + geom_cross()

quickcor(mtcars, cor.test = TRUE) + geom_star(n = 5)

quickcor(mtcars, cor.test = TRUE) + geom_colour() + geom_number(aes(num = r))

quickcor(mtcars, cor.test = TRUE) +
  geom_square(data = get_data(type = "lower", show.diag = FALSE)) +
  geom_mark(data = get_data(type = "upper", show.diag = FALSE)) +
  geom_abline(slope = -1, intercept = 12)

ggcor()函数

ggcor()函数不支持内部计算,所有首先要通过fortify_cor()或者as_cor_tbl()函数将数据转换成需要的格式。高能预警:默认出图可能有点丑,不想自己动手处理的可以跳过了。

args(ggcor)
## function (data, mapping = NULL, axis.x.position = "auto", axis.y.position = "auto", 
##     axis.label.drop = TRUE) 
## NULL
df01 <- fortify_cor(mtcars)
ggcor(df01) # 最基本的情况

ggcor(df01) + geom_square()

ggcor(df01) + 
  add_grid() +
  geom_square(aes(r0 = r, fill = r)) +
  scale_fill_gradient2n() +
  coord_fixed(expand = FALSE) +
  theme_cor() ## 这个例子还原了quickcor()在内部做了哪些工作。

ggcor(df01) + geom_point(aes(size = r)) +
  scale_size_area(max_size = 10)

样本分组的可视化

从数据结构来说,样本分组就是根据某个指示变量,将原数据表按行拆成子表,然后对子表做相关性分析。我们用iris数据集来看看具体怎么回事。

iris数据集的最后一列(Species)是样本分类信息,分别是三种不同的鸾尾花,当然我们可以直接把这个数据集当成一个表,来分析几个变量之间的相关性。

DT::datatable(iris)
quickcor(iris[-5]) + geom_colour()

但是,对于有些研究来说,可能想知道每一种鸾尾花的几个测量指标直接的相关性,那么对于大多数的同类的包,只能分开处理,然后手动拼图,我在ggcor中整合了这种情况。group参数必须是一个和原数据等长的向量。

从结果来看,样本分组之后相关性和混在一起的情况完全不一样,所以自我安慰下这个是有用的。

quickcor(iris[-5], group = iris[[5]]) + geom_colour() +
  facet_wrap(~group)

quickcor(iris[-5], type = "upper", show.diag = TRUE, group = iris[[5]]) + 
  geom_colour() +
  facet_wrap(~group)

grp <- mtcars$gear
quickcor(mtcars[-10], type = "lower", group = grp) + 
  geom_colour() +
  add_diaglab(hjust = 0) +
  expand_axis(x = 12) +
  remove_axis("x") +
  facet_wrap(~group)

以为只能对对数的相关系数矩阵分面?那是不可能的,我们再来看一个虚拟的例子。

group <- rep(paste0("sample", 1:2), 12)
quickcor(varechem[ , 1:6], varespec[ , 1:12],  group = group) + 
  geom_colour() +
  facet_grid(rows = vars(group))

变化(六)

对于mantel test的数据处理变化比较大,原来函数设计是以画图为基础,现在把数据处理进一步泛化,应用范围更广。

fortify_mantel()

这里要说明两个情况:第一,新版本物种数据(spec)和环境数据(env)必须是data.frame或者其它可以转换成data.frame的对象,原版本对于list的支持已近放弃了;第二,当不设置其它参数的时候,是把spec当成一个物种(一类),然后计算和env中每一列的相关性。

args(fortify_mantel)
## function (spec, env, group = NULL, env.ctrl = NULL, mantel.fun = "mantel", 
##     ...) 
## NULL
fortify_mantel(varespec, varechem)
## # A tibble: 14 x 4
##    spec  env             r p.value
##  * <chr> <chr>       <dbl>   <dbl>
##  1 spec  N         0.152     0.036
##  2 spec  P         0.221     0.009
##  3 spec  K         0.213     0.009
##  4 spec  Ca        0.0961    0.151
##  5 spec  Mg        0.0634    0.238
##  6 spec  S         0.137     0.062
##  7 spec  Al        0.278     0.001
##  8 spec  Fe        0.167     0.049
##  9 spec  Mn        0.221     0.006
## 10 spec  Zn        0.0400    0.31 
## 11 spec  Mo        0.00330   0.463
## 12 spec  Baresoil  0.145     0.044
## 13 spec  Humdepth  0.188     0.016
## 14 spec  pH       -0.0282    0.602

列分组

这里主要是两个参数:spec.seletenv.select。前者对spec数据框进行列分组,后者对env数据框进行列分组。这两个参数需要列表(list)形式,最好对每一组提供一个有意义的名字。有很多人问我,这个列表里面的数字是什么?其实就是你要把spec的那几列归成一组,数字就是列的序号。(若是你看不懂1:5是什么,那么我告诉你,你需要找一本R入门的书看看。

## 若spec的1-5列是spec01, 6-12列是spec02
fortify_mantel(varespec, varechem, 
               spec.select = list(spec01 = 1:5, spec02 = 6:12))
## # A tibble: 28 x 4
##    spec   env          r p.value
##  * <chr>  <chr>    <dbl>   <dbl>
##  1 spec01 N      0.256     0.025
##  2 spec01 P      0.137     0.086
##  3 spec01 K      0.399     0.004
##  4 spec01 Ca     0.00430   0.456
##  5 spec01 Mg     0.0254    0.36 
##  6 spec01 S      0.282     0.013
##  7 spec01 Al     0.313     0.009
##  8 spec01 Fe     0.0806    0.276
##  9 spec01 Mn    -0.0575    0.63 
## 10 spec01 Zn    -0.119     0.798
## # … with 18 more rows
## 若spec的1-5列是spec01, 6-12列是spec02
## 若env的1-5是env01,6-10是env02,11-14是env03
fortify_mantel(varespec, varechem, 
               spec.select = list(spec01 = 1:5, spec02 = 6:12),
               env.select = list(env01 = 1:5, env02 = 6:10, env03 = 11:14))
## # A tibble: 6 x 4
##   spec   env         r p.value
## * <chr>  <chr>   <dbl>   <dbl>
## 1 spec01 env01  0.0556   0.289
## 2 spec01 env02  0.285    0.017
## 3 spec01 env03 -0.0820   0.752
## 4 spec02 env01  0.0371   0.309
## 5 spec02 env02  0.0928   0.107
## 6 spec02 env03  0.0301   0.302
## 若spec和env都是单列一组
##++++++++++++++++++++++++++++++++++++++++++##
## 这里有一个很常见的报错,spec默认bray距离,不能出现行的和为0,否则会报错。
## 所以,这里设置spec.dist.method为欧几里得距离算法。
## 当然,为了计算速度更快,mantel.fun使用了ade4包的mantel.randtest()函数
##++++++++++++++++++++++++++++++++++++++++++##
sp.select <- as.list(setNames(names(varespec), names(varespec)))
fortify_mantel(varespec, varechem, spec.select = sp.select,
               spec.dist.method = "euclidean", mantel.fun = "mantel.randtest")
## # A tibble: 616 x 4
##    spec     env         r p.value
##  * <chr>    <chr>   <dbl>   <dbl>
##  1 Callvulg N      0.0129   0.397
##  2 Callvulg P     -0.129    0.969
##  3 Callvulg K      0.418    0.03 
##  4 Callvulg Ca    -0.126    0.887
##  5 Callvulg Mg    -0.0283   0.443
##  6 Callvulg S      0.299    0.047
##  7 Callvulg Al     0.199    0.064
##  8 Callvulg Fe     0.0713   0.168
##  9 Callvulg Mn    -0.0607   0.472
## 10 Callvulg Zn    -0.113    0.842
## # … with 606 more rows

样本分组

和相关性分析类似,有时候在不同的时间,不同的地点采样,就需要把样本分成几组分别分析。整个处理方法就是要指定样本分类的标签(group参数)。数据限制,以下例子均是虚构,完全随机分类。

set.seed(20191224)
sam_grp <- sample(paste0("sample", 1:3), 24, replace = TRUE)
fortify_mantel(varespec, varechem, group = sam_grp)
## # A tibble: 42 x 5
##    spec  env         r p.value group  
##    <chr> <chr>   <dbl>   <dbl> <chr>  
##  1 spec  N      0.0917   0.306 sample1
##  2 spec  P      0.0853   0.253 sample1
##  3 spec  K      0.205    0.133 sample1
##  4 spec  Ca    -0.160    0.759 sample1
##  5 spec  Mg    -0.0812   0.639 sample1
##  6 spec  S     -0.138    0.745 sample1
##  7 spec  Al     0.561    0.006 sample1
##  8 spec  Fe     0.463    0.02  sample1
##  9 spec  Mn     0.218    0.155 sample1
## 10 spec  Zn    -0.0325   0.538 sample1
## # … with 32 more rows
## 同时列分组和样本分组
fortify_mantel(varespec, varechem, group = sam_grp,
               spec.select = list(spec01 = 1:5, spec02 = 6:12))
## # A tibble: 84 x 5
##    spec   env         r p.value group  
##    <chr>  <chr>   <dbl>   <dbl> <chr>  
##  1 spec01 N      0.312    0.093 sample1
##  2 spec01 P      0.317    0.054 sample1
##  3 spec01 K      0.400    0.062 sample1
##  4 spec01 Ca     0.178    0.235 sample1
##  5 spec01 Mg     0.0880   0.282 sample1
##  6 spec01 S      0.0795   0.36  sample1
##  7 spec01 Al     0.308    0.121 sample1
##  8 spec01 Fe    -0.0195   0.407 sample1
##  9 spec01 Mn     0.0771   0.409 sample1
## 10 spec01 Zn    -0.0730   0.535 sample1
## # … with 74 more rows

mantel test 可视化

最简单的,把mantel test结果当成相关系数来处理,可以直接把数据结果传递给quickcor()函数。

mantel01 <- fortify_mantel(varespec, varechem, 
                           spec.select = list(spec01 = 1:5, 
                                              spec02 = 6:12,
                                              spec03 = 7:18,
                                              spec04 = 20:29,
                                              spec05 = 30:44))
quickcor(mantel01, legend.title = "Mantel's r") + 
  geom_square() + geom_cross()

mantel02 <- fortify_mantel(varespec, varechem, group = sam_grp,
                           spec.select = list(spec01 = 1:5, 
                                              spec02 = 6:12,
                                              spec03 = 7:18,
                                              spec04 = 20:29,
                                              spec05 = 30:44),
                           mantel.fun = "mantel.randtest")
quickcor(mantel02, legend.title = "Mantel's r") + 
  geom_colour() + geom_cross() + facet_grid(rows = vars(group))

science 组合图

作为ggcor最大的不一样,组合图的部分还是需要重点介绍。函数还是add_link(),不过我已经放弃在函数内部猜用户的心思了,每个人的需求不一样,做得越多反而越麻烦,所以新版本的add_link()函数仅仅把线条画出来,至于你是需要什么线型、什么颜色等等,都需要自己去整。

## 默认出图
mantel03 <- fortify_mantel(varespec, varechem, 
                         spec.select = list(1:10, 5:14, 7:22, 9:32))
quickcor(varechem, type = "upper") + geom_square() + 
  add_link(mantel03, diag.label = TRUE) +
  add_diaglab() + remove_axis("x")

## 调整映射
library(dplyr)
mantel031 <- mantel03 %>% 
  mutate(r = cut(r, breaks = c(-Inf, 0.25, 0.5, Inf), 
                 labels = c("<0.25", "0.25-0.5", ">=0.5"),
                 right = FALSE),
         p.value = cut(p.value, breaks = c(-Inf, 0.001, 0.01, 0.05, Inf),
                       labels = c("<0.001", "0.001-0.01", "0.01-0.05", ">=0.05"),
                       right = FALSE))
quickcor(varechem, type = "upper") + geom_square() + 
  add_link(mantel031, mapping = aes(colour = p.value, size = r),
           diag.label = TRUE) +
  scale_size_manual(values = c(0.5, 1.5, 3)) +
  add_diaglab() + remove_axis("x")

## 剔除不显著的
mantel032 <- mantel03 %>% filter(p.value < 0.05) %>% 
  mutate(r = cut(r, breaks = c(-Inf, 0.25, 0.5, Inf), 
                 labels = c("<0.25", "0.25-0.5", ">=0.5"),
                 right = FALSE),
         p.value = cut(p.value, breaks = c(-Inf, 0.001, 0.01, 0.05, Inf),
                       labels = c("<0.001", "0.001-0.01", "0.01-0.05", ">=0.05"),
                       right = FALSE))
quickcor(varechem, type = "upper") + geom_square() + 
  add_link(mantel032, mapping = aes(colour = p.value, size = r),
           diag.label = TRUE) +
  scale_size_manual(values = c(0.5, 1.5, 3)) +
  add_diaglab() + remove_axis("x")

其它数据

df02 <- fortify_cor(varespec[1:4], varechem, cor.test = TRUE) %>% 
  mutate(p.value = cut(p.value, breaks = c(-Inf, 0.05, Inf),
                       labels = c("<0.05", ">=0.05"), right = FALSE))
quickcor(varechem, type = "upper") + geom_square() + 
  add_link(df02, mapping = aes(colour = p.value, size = r), alpha = 0.5,
           spec.key = "idy", env.key = "idx",
           diag.label = TRUE) +
  scale_size_continuous(range = c(0.2, 1.5))  +
  scale_color_manual(values = c("#D95F02", "#CCCCCC")) +
  add_diaglab() + 
  guides(colour = guide_legend(title = "p value", order = 2), 
         size = guide_legend(title = "correlation", order = 3)) +
  remove_axis("x")