厚缊

诹图——ggcor简介(六)

厚缊 / 2020-02-06


Mantel 检验是相关性检验的一种泛化,本质上是比较两个矩阵的相关性。ggcor只是对已有包(veganade4)相关函数做了再封装,没有提供这方面的算法实现。

基础函数

fortify_mantel()mantel_test()函数是ggcor包封装的两个mantel 检验函数,其中mantel_test()是最基本的函数,fortify_mantel()依赖于mantel_test()函数。

mantel 检验本质上是对称的,并不会严格区分物种(special)、环境(environment)和控制环境(control environment),但是为了更好地区别这两者,抑或说是一种设计惯性,让我在mantel_test()函数中对参数命名做了区分。

必须的两个参数是specenv,要求必须是data frame对象或者其它可以转成data frame对象的结构。mantel.fun是指定检验函数,目前支持四种,默认是vegan::mantel(),详细的请看帮助文档。

需要注意的是,默认情况下spec使用“bray”距离算法,env使用“euclidean”距离算法,若要自行指定可分别通过spec.dist.methodenv.dist.method参数设置。

library(ggcor)
args(mantel_test)
## function (spec, env, env.ctrl = NULL, mantel.fun = "mantel", 
##     spec.select = NULL, env.select = NULL, spec.dist.method = "bray", 
##     env.dist.method = "euclidean", ...) 
## NULL

基本用法

除了数据,不设置任何参数时是把spec的所有列当成一个群落,然后分别和env的每一列做mantel 检验。

data("varespec", package = "vegan")
data("varechem", package = "vegan")
mantel_test(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.012
##  4 spec  Ca        0.0961    0.137
##  5 spec  Mg        0.0634    0.243
##  6 spec  S         0.137     0.053
##  7 spec  Al        0.278     0.003
##  8 spec  Fe        0.167     0.046
##  9 spec  Mn        0.221     0.007
## 10 spec  Zn        0.0400    0.314
## 11 spec  Mo        0.00330   0.453
## 12 spec  Baresoil  0.145     0.05 
## 13 spec  Humdepth  0.188     0.025
## 14 spec  pH       -0.0282    0.556

上面的例子等价于:

mantel_test(varespec, varechem, 
            spec.select = list(spec = 1:length(varespec)))
## # A tibble: 14 x 4
##    spec  env             r p.value
##  * <chr> <chr>       <dbl>   <dbl>
##  1 spec  N         0.152     0.033
##  2 spec  P         0.221     0.004
##  3 spec  K         0.213     0.007
##  4 spec  Ca        0.0961    0.125
##  5 spec  Mg        0.0634    0.257
##  6 spec  S         0.137     0.047
##  7 spec  Al        0.278     0.003
##  8 spec  Fe        0.167     0.049
##  9 spec  Mn        0.221     0.018
## 10 spec  Zn        0.0400    0.297
## 11 spec  Mo        0.00330   0.456
## 12 spec  Baresoil  0.145     0.045
## 13 spec  Humdepth  0.188     0.019
## 14 spec  pH       -0.0282    0.598

多个群落

多个群落的情况可能更常见,可以通过spec.select参数来控制。这里可以大概说说spec.select的作用机制,该参数指定了spec中哪些列构成一个群落,下面的例子中spec01 = 1:5表示spec的1-5列构成spec01群落,其它的以此类推。除了默认情况,spec可以有冗余变量,因为没有没有spec.select选中的列自动丢弃。

mantel_test(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.02 
##  2 spec01 P      0.137     0.084
##  3 spec01 K      0.399     0.004
##  4 spec01 Ca     0.00430   0.416
##  5 spec01 Mg     0.0254    0.354
##  6 spec01 S      0.282     0.025
##  7 spec01 Al     0.313     0.014
##  8 spec01 Fe     0.0806    0.268
##  9 spec01 Mn    -0.0575    0.632
## 10 spec01 Zn    -0.119     0.796
## # … with 18 more rows

除了可以通过列序号指定外,还可以通过列名来指定。

name <- names(varespec)
name
##  [1] "Callvulg" "Empenigr" "Rhodtome" "Vaccmyrt" "Vaccviti" "Pinusylv"
##  [7] "Descflex" "Betupube" "Vacculig" "Diphcomp" "Dicrsp"   "Dicrfusc"
## [13] "Dicrpoly" "Hylosple" "Pleuschr" "Polypili" "Polyjuni" "Polycomm"
## [19] "Pohlnuta" "Ptilcili" "Barbhatc" "Cladarbu" "Cladrang" "Cladstel"
## [25] "Cladunci" "Cladcocc" "Cladcorn" "Cladgrac" "Cladfimb" "Cladcris"
## [31] "Cladchlo" "Cladbotr" "Cladamau" "Cladsp"   "Cetreric" "Cetrisla"
## [37] "Flavniva" "Nepharct" "Stersp"   "Peltapht" "Icmaeric" "Cladcerv"
## [43] "Claddefo" "Cladphyl"
mantel_test(varespec, varechem, 
            spec.select = list(spec01 = name[1:5],
                               spec02 = name[6:10]))
## # A tibble: 28 x 4
##    spec   env          r p.value
##  * <chr>  <chr>    <dbl>   <dbl>
##  1 spec01 N      0.256     0.015
##  2 spec01 P      0.137     0.085
##  3 spec01 K      0.399     0.002
##  4 spec01 Ca     0.00430   0.461
##  5 spec01 Mg     0.0254    0.344
##  6 spec01 S      0.282     0.014
##  7 spec01 Al     0.313     0.01 
##  8 spec01 Fe     0.0806    0.273
##  9 spec01 Mn    -0.0575    0.583
## 10 spec01 Zn    -0.119     0.775
## # … with 18 more rows

环境变量分类

spec可以由多列构成一个群落,env就可以由多列构成一类,设置方式、原理均与spec类似。

mantel_test(varespec, varechem, 
            env.select = list(env01 = 1:7,
                              env02 = 8:14))
## # A tibble: 2 x 4
##   spec  env       r p.value
## * <chr> <chr> <dbl>   <dbl>
## 1 spec  env01 0.174   0.029
## 2 spec  env02 0.241   0.008

specenv是可以同时分类的。

mantel_test(varespec, varechem, 
            spec.select = list(spec01 = 1:5,
                               spec02 = 6:12),
            env.select = list(env01 = 1:7,
                              env02 = 8:14))
## # A tibble: 4 x 4
##   spec   env        r p.value
## * <chr>  <chr>  <dbl>   <dbl>
## 1 spec01 env01 0.163    0.103
## 2 spec01 env02 0.0416   0.379
## 3 spec02 env01 0.0633   0.214
## 4 spec02 env02 0.154    0.036

偏mantel 检验

目前来说,偏mantel 检验的实现还是比较弱,,只能设置一个统一的控制环境,不能给每一对specenv组合指定不同的控制环境(env.ctrl)。

mantel_test(varespec, varechem, env.ctrl = varechem[7:14],
            spec.select = list(spec01 = 1:5,
                               spec02 = 6:12),
            env.select = list(env01 = 1:3,
                              env02 = 4:6),
            mantel.fun = "mantel.partial")
## # A tibble: 4 x 4
##   spec   env          r p.value
## * <chr>  <chr>    <dbl>   <dbl>
## 1 spec01 env01  0.378     0.004
## 2 spec01 env02 -0.00614   0.493
## 3 spec02 env01 -0.110     0.9  
## 4 spec02 env02  0.0457    0.271

其他参数

当数据量比较大的时候,建议把mantel.fun参数设置为“mantel.randtest”,这时调用的是ade4包中的mantel.randtest()进行mantel 检验,底层代码是C语言,计算速度比较快。

library(microbenchmark)
microbenchmark(
  mantel1 = mantel_test(varespec, varechem, 
                        mantel.fun = "mantel"),
  mantel2 = mantel_test(varespec, varechem, mantel.fun = "mantel",
              parallel = 2),
  mantel.rand = mantel_test(varespec, varechem, 
                            mantel.fun = "mantel.randtest"),
  mantel.r = mantel_test(varespec, varechem, 
                         mantel.fun = "mantel.rtest"),
  times = 10
)
## Unit: milliseconds
##         expr        min         lq       mean     median         uq
##      mantel1  745.68224  747.98636  800.65029  778.13250  795.17616
##      mantel2 1282.25915 1372.56686 1428.22702 1414.87283 1477.89454
##  mantel.rand   59.89069   68.02138   74.88553   70.07595   71.36375
##     mantel.r  190.67209  192.19817  203.64338  200.74595  214.00982
##        max neval
##   940.4380    10
##  1626.1158    10
##   129.1531    10
##   227.8423    10

其他mantel 检验的函数细节控制需要根据不同的mantel.fun进行针对性设置。下面的例子说明了如何控制置换检验的次数。

args(vegan::mantel) 
## function (xdis, ydis, method = "pearson", permutations = 999, 
##     strata = NULL, na.rm = FALSE, parallel = getOption("mc.cores")) 
## NULL
args(ade4::mantel.randtest)
## function (m1, m2, nrepet = 999, ...) 
## NULL
mantel_test(varespec, varechem, mantel.fun = "mantel",
            permutations = 2000) ## 默认是999次
## # A tibble: 14 x 4
##    spec  env             r p.value
##  * <chr> <chr>       <dbl>   <dbl>
##  1 spec  N         0.152   0.0390 
##  2 spec  P         0.221   0.00750
##  3 spec  K         0.213   0.0110 
##  4 spec  Ca        0.0961  0.130  
##  5 spec  Mg        0.0634  0.228  
##  6 spec  S         0.137   0.0535 
##  7 spec  Al        0.278   0.00200
##  8 spec  Fe        0.167   0.0485 
##  9 spec  Mn        0.221   0.0115 
## 10 spec  Zn        0.0400  0.308  
## 11 spec  Mo        0.00330 0.484  
## 12 spec  Baresoil  0.145   0.0485 
## 13 spec  Humdepth  0.188   0.0215 
## 14 spec  pH       -0.0282  0.606
mantel_test(varespec, varechem, mantel.fun = "mantel.randtest",
            nrepet = 2000) ## 默认是999次
## # A tibble: 14 x 4
##    spec  env             r p.value
##  * <chr> <chr>       <dbl>   <dbl>
##  1 spec  N         0.152   0.0220 
##  2 spec  P         0.221   0.00350
##  3 spec  K         0.213   0.0105 
##  4 spec  Ca        0.0961  0.136  
##  5 spec  Mg        0.0634  0.236  
##  6 spec  S         0.137   0.0545 
##  7 spec  Al        0.278   0.00200
##  8 spec  Fe        0.167   0.0435 
##  9 spec  Mn        0.221   0.00850
## 10 spec  Zn        0.0400  0.328  
## 11 spec  Mo        0.00330 0.454  
## 12 spec  Baresoil  0.145   0.0355 
## 13 spec  Humdepth  0.188   0.0220 
## 14 spec  pH       -0.0282  0.598

fortify_mantel()函数

fortify_mantel()函数是mantel_test()函数的再包装,主要目的是为了增加样本分组的情况,也就是针对不同的样本组分别做mantel 检验。

set.seed(20200206)
sample_group <- sample(LETTERS[1:2], 24, replace = TRUE)
fortify_mantel(varespec, varechem, group = sample_group)
## # A tibble: 28 x 5
##    spec  env          r p.value .group
##    <chr> <chr>    <dbl>   <dbl> <chr> 
##  1 spec  N     -0.0237    0.531 A     
##  2 spec  P     -0.0369    0.545 A     
##  3 spec  K      0.0828    0.246 A     
##  4 spec  Ca    -0.266     0.973 A     
##  5 spec  Mg    -0.174     0.893 A     
##  6 spec  S     -0.00476   0.466 A     
##  7 spec  Al     0.541     0.001 A     
##  8 spec  Fe     0.327     0.023 A     
##  9 spec  Mn     0.202     0.099 A     
## 10 spec  Zn    -0.0600    0.637 A     
## # … with 18 more rows

为了帮助大家更好的理解样本分组,我们来一步步的看看fortify_mantel()是如何做到的。

  1. 根据group参数(注意,必须是和specenv行数相等的向量)将specenv拆分成列表;

  2. 分别对specenv列表里面的对应元素(spec列表的第一个元素对应env列表的第一个元素,以此类推)调用mantel_test()函数进行mantel 检验。

  3. 将每一组检验的结果按行堆叠起来,并增加.group列。

从这里可以看出,当不设置group参数(或者说不进行样本分组)时,fortify_mantel()完全等价于mantel_test()函数,当设置了group参数时,fortify_mantel()的其它参数(除了env.ctrl,等下另说)的均和mantel_test()相同。

mantel_test(varespec, varechem)
## # A tibble: 14 x 4
##    spec  env             r p.value
##  * <chr> <chr>       <dbl>   <dbl>
##  1 spec  N         0.152     0.033
##  2 spec  P         0.221     0.004
##  3 spec  K         0.213     0.02 
##  4 spec  Ca        0.0961    0.128
##  5 spec  Mg        0.0634    0.232
##  6 spec  S         0.137     0.06 
##  7 spec  Al        0.278     0.004
##  8 spec  Fe        0.167     0.044
##  9 spec  Mn        0.221     0.013
## 10 spec  Zn        0.0400    0.327
## 11 spec  Mo        0.00330   0.488
## 12 spec  Baresoil  0.145     0.043
## 13 spec  Humdepth  0.188     0.028
## 14 spec  pH       -0.0282    0.612
fortify_mantel(varespec, varechem)
## # A tibble: 14 x 4
##    spec  env             r p.value
##  * <chr> <chr>       <dbl>   <dbl>
##  1 spec  N         0.152     0.03 
##  2 spec  P         0.221     0.009
##  3 spec  K         0.213     0.009
##  4 spec  Ca        0.0961    0.124
##  5 spec  Mg        0.0634    0.239
##  6 spec  S         0.137     0.068
##  7 spec  Al        0.278     0.002
##  8 spec  Fe        0.167     0.05 
##  9 spec  Mn        0.221     0.017
## 10 spec  Zn        0.0400    0.301
## 11 spec  Mo        0.00330   0.481
## 12 spec  Baresoil  0.145     0.047
## 13 spec  Humdepth  0.188     0.016
## 14 spec  pH       -0.0282    0.599

带控制环境和样本分组的mantel 检验

当需要进行样本分组时,每个分组的样本数量可能并不是相同的,这个时候至少需要给每个样本组一个控制环境,fortify_mantel()是支持的。

table(sample_group)
## sample_group
##  A  B 
## 11 13
fortify_mantel(varespec, varechem, 
               env.ctrl = list(A = varechem[1:11, 1:3],
                               B = varechem[12:24, 12:14]),
               group = sample_group, mantel.fun = "mantel.partial")
## # A tibble: 28 x 5
##    spec  env         r p.value .group
##    <chr> <chr>   <dbl>   <dbl> <chr> 
##  1 spec  N     -0.0319   0.545 A     
##  2 spec  P     -0.0518   0.618 A     
##  3 spec  K      0.0586   0.295 A     
##  4 spec  Ca    -0.148    0.84  A     
##  5 spec  Mg    -0.0819   0.685 A     
##  6 spec  S     -0.0113   0.476 A     
##  7 spec  Al     0.546    0.001 A     
##  8 spec  Fe     0.297    0.03  A     
##  9 spec  Mn     0.139    0.184 A     
## 10 spec  Zn    -0.114    0.783 A     
## # … with 18 more rows

把mantel 检验结果转成cor_tbl

as_cor_tbl()直接支持这种转换。

mantel_test(varespec, varechem, 
            spec.select = list(spec01 = 1:5,
                               spec02 = 6:12,
                               spec03 = 13:17,
                               spec04 = 18:28)) %>% 
  as_cor_tbl() %>% 
  quickcor() + geom_colour()

mantel_test(varespec, varechem, 
            spec.select = list(spec01 = 1:5,
                               spec02 = 6:12,
                               spec03 = 13:17,
                               spec04 = 18:28)) %>% 
  as_cor_tbl(byrow = FALSE) %>% 
  quickcor() + geom_colour()

结果的格式化输出

因为mantel 检验结果和cor_tbl之间是可以互相转化的,所以display_cor()也同样支持mantel 检验结果。

mantel_test(varespec, varechem, 
            spec.select = list(spec01 = 1:5,
                               spec02 = 6:12,
                               spec03 = 13:17,
                               spec04 = 18:28)) %>% 
  display_cor(byrow = FALSE) %>% 
  DT::datatable()

小结

目前来说,限于没有专业背景,ggcor对于mantel 检验就做了这么多工作,若有更好的建议,可以直接联系我。

下期预告:science 组合图。