这一期会很简单,主要是提供一种mantel检验的替代方案——procrutes 检验。procrutes 检验在图像识别领域用的比较多,用来比较两组数据之间的相似性,本质上说是通过平移、旋转、缩放等一系列变换,使原数据集的点匹配到目标数据集上,并使其和对应点的误差平方和最小。和mantel类似,p值也是通过置换检验得到的。原理的东西确实懂得不多,只能班门弄斧,强行解释。
这部分内容应该是0.9.4就加进去了,一直没有介绍,若是0.9.4版本,可以不用更新。
if(!require(ggcor)) {
devtools::install_github('houyunhuang/ggcor')
}
## 若之前安装过老版本,加上`force = TRUE`参数
packageVersion('ggcor')
## [1] '0.9.5'
ggcor
并不会像vegan
、ade4
包一样,详细的分析一组procrutes 检验的结果,而是用来处理多组之间比较和检验。默认情况下,spec
整体当成一组,env
是每列一组,然后分析每个spec
组和每个env
组之间的关系,也就是说默认情况下是比较的spec
整体和env
每个单列之间的相似性。procrutes_test()
和mantel_test()
函数之所以有如此变态的设置,是因为science的那个组合图,但是这样也会误导,让人觉得这两种检验就该这样使用。
library(ggcor)
args(procrutes_test)
## function (spec, env, procrutes.fun = 'protest', spec.select = NULL,
## env.select = NULL, spec.pre.fun = 'identity',
## spec.pre.params = list(),
## env.pre.fun = spec.pre.fun,
## env.pre.params = spec.pre.params,
## ...)
## NULL
data('varespec', package = 'vegan')
data('varechem', package = 'vegan')
procrutes_test(varespec, varechem)
## # A tibble: 14 x 4
## spec env r p.value
## * <chr> <chr> <dbl> <dbl>
## 1 spec N 0.326 0.072
## 2 spec P 0.329 0.041
## 3 spec K 0.301 0.092
## 4 spec Ca 0.324 0.059
## 5 spec Mg 0.302 0.087
## 6 spec S 0.251 0.21
## 7 spec Al 0.452 0.003
## 8 spec Fe 0.410 0.01
## 9 spec Mn 0.443 0.003
## 10 spec Zn 0.269 0.152
## 11 spec Mo 0.139 0.76
## 12 spec Baresoil 0.428 0.006
## 13 spec Humdepth 0.436 0.008
## 14 spec pH 0.358 0.016
若我们想自己设置spec
中的分组顺序,这和mantel_test()
完全一致,通过spec.select
参数设置,同样的道理也适用于env
。再次强调,每个分组都要指定名字,否则内部直接自动命名。
## 自定义spec分组
procrutes_test(varespec, varechem,
spec.select = list(A = 1:7, B = 8:40))
## # A tibble: 28 x 4
## spec env r p.value
## * <chr> <chr> <dbl> <dbl>
## 1 A N 0.301 0.073
## 2 A P 0.107 0.926
## 3 A K 0.260 0.191
## 4 A Ca 0.219 0.314
## 5 A Mg 0.205 0.378
## 6 A S 0.211 0.415
## 7 A Al 0.392 0.006
## 8 A Fe 0.309 0.08
## 9 A Mn 0.202 0.424
## 10 A Zn 0.104 0.895
## # … with 18 more rows
## 自定义spec和env分组
procrutes_test(varespec, varechem,
spec.select = list(A = 1:7, B = 8:40),
env.select = list(Env01 = 1:7, Env02 = 8:14))
## # A tibble: 4 x 4
## spec env r p.value
## * <chr> <chr> <dbl> <dbl>
## 1 A Env01 0.382 0.063
## 2 A Env02 0.367 0.119
## 3 B Env01 0.476 0.007
## 4 B Env02 0.524 0.001
procrutes检验比mantel检验复杂的地方是,mantel是先进行距离变换,而procrutes可以做任何变换,只要保证输出结果是矩阵即可。目前我看到的有主成份降维和monoMDS处理。procrutes_test()
函数提供了非常灵活的接口,支持自定义预处理方式。默认情况下,spec
和env
的预处理函数完全一样,都是identity
,即不做任何变换。
我们先来看看主成份变换的例子。dudi_pca()
函数是ggcor
封装的ade4::dudi.pca()
函数,若是要自定义主成份分析参数,可以通过spec.pre.params
处理。
procrutes_test(varespec, varechem,
spec.select = list(A = 1:7, B = 8:40),
env.select = list(Env01 = 1:7, Env02 = 8:14),
spec.pre.fun = 'dudi_pca')
## # A tibble: 4 x 4
## spec env r p.value
## * <chr> <chr> <dbl> <dbl>
## 1 A Env01 0.454 0.106
## 2 A Env02 0.445 0.252
## 3 B Env01 0.496 0.036
## 4 B Env02 0.557 0.003
## 求pca之前先中心标准化
procrutes_test(varespec, varechem,
spec.select = list(A = 1:7, B = 8:40),
env.select = list(Env01 = 1:7, Env02 = 8:14),
spec.pre.fun = 'dudi_pca',
spec.pre.params = list(center = TRUE))
## # A tibble: 4 x 4
## spec env r p.value
## * <chr> <chr> <dbl> <dbl>
## 1 A Env01 0.454 0.096
## 2 A Env02 0.445 0.245
## 3 B Env01 0.496 0.044
## 4 B Env02 0.557 0.001
monoMDS
变换也很简单,只需要调整下spec.pre.fun
参数就行。
procrutes_test(varespec, varechem,
spec.select = list(A = 1:7, B = 8:40),
env.select = list(Env01 = 1:7, Env02 = 8:14),
spec.pre.fun = 'mono_mds')
## # A tibble: 4 x 4
## spec env r p.value
## * <chr> <chr> <dbl> <dbl>
## 1 A Env01 0.416 0.029
## 2 A Env02 0.358 0.091
## 3 B Env01 0.521 0.006
## 4 B Env02 0.532 0.002
也有可能,ggcor
直接支持的上述三种预处理方式都不能满足你的需求,那么你需要自定义预处理函数。这里我们定义一个行极值标准化的预处理函数作为例子,这个只是表面可以这样定义,不具备理论意义。
row_scale <- function(x, na.rm = TRUE) {
if(!is.matrix(x))
x <- as.matrix(x)
row.max <- apply(x, 1, max, na.rm = na.rm)
row.min <- apply(x, 1, min, na.rm = na.rm)
(x - row.min) / (row.max - row.min)
}
procrutes_test(varespec, varechem,
spec.select = list(A = 1:7, B = 8:40),
env.select = list(Env01 = 1:7, Env02 = 8:14),
spec.pre.fun = 'row_scale')
## # A tibble: 4 x 4
## spec env r p.value
## * <chr> <chr> <dbl> <dbl>
## 1 A Env01 0.359 0.109
## 2 A Env02 0.320 0.389
## 3 B Env01 0.474 0.006
## 4 B Env02 0.512 0.003
既然procrutes检验是作为mantel检验的替代,那么当然也是可以很方便的来做组合图的,看个例子,不啰嗦了。
library(ggplot2)
corr <- fortify_cor(varechem, type = 'upper') ## 只能是上下三角 我们保留上三角
procrutes <- procrutes_test(varespec, varechem,
spec.select = list(Spec01 = 1:7,
Spec02 = 8:18,
Spec03 = 19:37,
Spec04 = 38:44),
spec.pre.fun = 'mono_mds') %>%
combination_layout(cor_tbl = corr) %>%
mutate(xend = xend + 1,
rd = cut(r, breaks = c(-Inf, 0.2, 0.4, Inf),
labels = c('< 0.2', '0.2 - 0.4', '>= 0.4')),
pd = cut(p.value, breaks = c(-Inf, 0.01, 0.05, Inf),
labels = c('< 0.01', '0.01 - 0.05', '>= 0.05')))
## 先关闭映射继承
options(ggcor.link.inherit.aes = FALSE)
quickcor(corr) + geom_square() +
geom_link(aes(colour = pd, size = rd), data = procrutes,
curvature = 0.05) +
geom_link_point(data = procrutes) +
geom_start_label(aes(x = x - 0.5), hjust = 1, data = procrutes) +
scale_size_manual(values = c(0.5, 1, 2)) +
scale_colour_manual(values = c('#D95F02', '#1B9E77', '#A2A2A288')) +
guides(size = guide_legend(title = 'Procrutes' r',
override.aes = list(colour = 'grey35'),
order = 2),
colour = guide_legend(title = 'Procrutes' p',
override.aes = list(size = 3),
order = 1),
fill = guide_colorbar(title = 'Pearson's r', order = 3)) +
expand_axis(x = -6)
procrutes检验的函数是完全凭感觉写出来的,符不符合使用习惯我也不清楚,若是有好的建议可以联系我。
联系客服