用向量做Mantel的几个问题

这几天有读者问我mental计算的几个问题，在此记录一下。

1. mantel test一般用距离矩阵来计算，vegan的mantel输入只能是距离矩阵。如果想用向量做mantel ，可以用ecodist包做，输入数据可以是向量的形式。
2. ecodist针对r<=0,r=0,r>=0分别输出了3个P值，不确定用哪个。我测试了一下发现r<=0时对应的P值和vegan中mantel结果的P值是一致的。因此可以用r<=0对应的P值，这也可以反推出vegan中mantel的原假设也是r<=0。 这一点在介绍ecodist的文章中已经添加： R——ecodist&MRM methods
3. 后来他又发现数据为435行时可以出结果，而当数据为704行的时候会报错： Error in mantel(A ~ B) : Matrix not square.

我测试了一下果然如此。把435随机换成其他几个数也会报错。 这时候开始有点意思了。难道435这个数存在什么特别之处么。 函数说明中没有提到这个报错，我在网上搜了一下也没有找到答案。

那就简单粗暴看函数源代码： https://github.com/cran/ecodist/blob/master/R/mantel.R

其中有这么几句：

1 m <- as.matrix(m)
2 n <- (1 + sqrt(1 + 8 * nrow(m)))/2  
3  if (abs(n - round(n)) > 1e-07) 
4    stop("Matrix not square.\n")

m为输入数据，通过一个公式得到n，将n和其整数部分进行比较，如果不相等则会报这个错。 接下来测试当输入的数据是多少行的时候不会报错：

1 s = c()
2  for (i in 1:1000){
3    n <- (1 + sqrt(1 + 8 * i))/2   
4    if (n - round(n) == 0 ){s = c(s,i)}
5  }
6s
7[1]   1   3   6  10  15  21  28  36  45  55  66  78  91 105 120 136 153 171 190
8[20] 210 231 253 276 300 325 351 378 406 435 465 496 528 561 595 630 666 703 741
9[39] 780 820 861 903 946 990

1000行以内只有44种情况不会报错，其中就有435。 结合这个结果和报错信息，我才突然发现原来输入数据的行数（1，3，6，10…）必须满足可以被转化为对称矩阵中上（或下）三角的形式才会计算结果。如435正好填满29*29的上（下）三角矩阵。其他数字得到的不是对称矩阵，因此会报错：Matrix not square。

所以ecodist用向量计算mantel还是有隐含的前提条件的。

如果数据不方便先转化为矩阵，那只能取特定的行数输入才能算mantel。

点分享

点点赞

点在看

一个环境工程专业却做生信分析的深井冰博士，深受拖延症的困扰。想给自己一点压力，争取能够不定期分享学到的生信小技能，亦或看文献过程中的一些笔记与小收获，记录生活中的杂七杂八。

目前能力有限，尚不能创造知识，只是知识的搬运工。