Mantel test的R实现

Mantel test 是对两个矩阵相关关系的检验，由Nathan Mantel在1976年提出。之所以抛开相关系数发展这样一种方法，是因为相关系数只能处理两列数据之间的相关性，而在面对两个矩阵之间的相关性时就束手无策。Mantel检验专治这种不服。

这种方法多用于生态学上，不同的样本case对应不同的变量，而不同的变量可以分属不同的类别，对case有不同角度的刻画。如基于不同植物种类数量可以建立样本间的两两距离矩阵，只需套用距离计算公式即可；不同样本的微生物clone序列，通过Unifrac方法也可以计算得到样本间距离矩阵；不同位置，两两间距离也可以用距离表示。所得到这些矩阵，如果希望验证两类描述间有没有相关关系，就非常有用了。比如我希望检验微生物群落是否和植被群落有对应关系，就可以将微生物Unifrac矩阵对植物的比如Bray-Curtis距离矩阵做个相关分析，由得到的结果得出自己的推论。这种方法的好处在于，不管你是什么数据，只要能计算有距离属性的值，都可以转化为距离矩阵进行分析。本人《空间方法在微生物生态学中的应用》（好像是生态学报2010年2月）一文中已有提及，在此延伸到实际应用。

Mantel test，顾名思义，是一种检验。既然是检验就得有原假设，它的原假设是两个矩阵见没有相关关系。检验过程如下：两个矩阵都对应展开，变量两列，计算相关系数（理论上什么相关系数都可以计算，但常用pearson相关系数），然后其中一列或两列同时置换，再计算一个值，permutation 成千上万次，看实际的r值在所得r值分布中的位置，如果跟随机置换得到的结果站队较近，则不大相关，如果远远比随机由此得到显著性。

Mantel测试测量通常包含距离测量的两个矩阵之间的相关性。 Mantel测试是一种测试空间自相关的方法。 在ade4库使用功能，我们可以在河进行Mantel检测要下载和加载这个库，输入install.packages（“ade4”），然后library(ade4)。 在其他R库中还有其他Mantel测试功能，我们对这个库的选择不应该被视为任何方式的认可。
让我们看一个例子。 我们的数据集中， 臭氧 ，含有在聚合一个月以上洛杉矶地区32的位置臭氧测量。 该数据集包括的站号（ 站 ），车站的纬度和经度（纬度和经度 ），最高8小时每日平均值（Av8top）的平均值。 我们将有兴趣测试是否距离较近的站相对于相距较远的站，臭氧测量的差异更小。 这些数据和其他空间数据集可以从伊利诺伊大学的空间分析实验室下载。 我们可以查看位置变量的摘要，以查看所考虑的位置范围。

ozone<-read.table("http://www.ats.ucla.edu/stat/r/faq/ozone.csv", sep=",", header=T)
head(ozone, n=10)

Station   Av8top      Lat       Lon
1       60 7.225806 34.13583 -117.9236
2       69 5.899194 34.17611 -118.3153
3       72 4.052885 33.82361 -118.1875
4       74 7.181452 34.19944 -118.5347
5       75 6.076613 34.06694 -117.7514
6       84 3.157258 33.92917 -118.2097
7       85 5.201613 34.01500 -118.0597
8       87 4.717742 34.06722 -118.2264
9       88 6.532258 34.08333 -118.1069
10      89 7.540323 34.38750 -118.5347

要运行Mantel测试，我们需要生成两个距离矩阵：一个包含空间距离，一个包含在给定点的测量结果之间的距离。 在空间距离矩阵中，靠近在一起的点对的条目低于远离的点对的条目。 在测量的结果矩阵中，具有类似结果的位置对的条目低于具有不同结果的点对的条目。 我们今天使用DIST功能做的。 Mantel测试函数将需要此“距离”类的对象。

station.dists <- dist(cbind(ozone$Lon, ozone$Lat))
ozone.dists <- dist(ozone$Av8top)

as.matrix(station.dists)[1:5, 1:5]

1         2         3         4         5
1 0.0000000 0.3937326 0.4088031 0.6144127 0.1854888
2 0.3937326 0.0000000 0.3749446 0.2206810 0.5743590
3 0.4088031 0.3749446 0.0000000 0.5116772 0.4994034
4 0.6144127 0.2206810 0.5116772 0.0000000 0.7944601
5 0.1854888 0.5743590 0.4994034 0.7944601 0.0000000

as.matrix(ozone.dists)[1:5, 1:5]

1        2        3        4        5
1 0.000000 1.326612 3.172921 0.044354 1.149193
2 1.326612 0.000000 1.846309 1.282258 0.177419
3 3.172921 1.846309 0.000000 3.128567 2.023728
4 0.044354 1.282258 3.128567 0.000000 1.104839
5 1.149193 0.177419 2.023728 1.104839 0.000000

这些是函数将要测试相关性的两个矩阵。 该测试包括计算矩阵中的条目的相关性，然后置换矩阵并且在每个置换下计算相同的测试统计量，并且将原始测试统计量与来自置换的测试统计量的分布进行比较以生成p值。 置换的数量定义了可以计算p值的精度。 进行Mantel检验的功能是mantel.rtest和必需的参数是两个距离矩阵。 排列数也可以由用户指定，但默认为99。

mantel.rtest(station.dists, ozone.dists, nrepet = 9999)

Monte-Carlo test
Observation: 0.1636308 
Call: mantel.rtest(m1 = station.dists, m2 = ozone.dists, nrepet = 9999)
Based on 9999 replicates
Simulated p-value: 0.0312 

基于这些结果，我们可以拒绝零假设，这两个矩阵，空间距离和臭氧距离是无关的 ****其中α= 0.05。 观察到的相关性，r = 0.1636308，表明矩阵条目是正相关的。 因此，在彼此靠近而彼此远离的站之间通常看到较小的臭氧差异。 注意，由于该测试是基于随机排列，所以相同的代码将总是达到相同的观察到的相关性，但很少具有相同的p值。

代码

# Mantel test -------------------------------------------------------------

# example from "ape"
library("ape", lib.loc="D:/Program Files/R/R-3.1.2/library")

q1 <- matrix(runif(36), nrow = 12)
q2 <- matrix(runif(36), nrow = 6)
mantel.test(q1, q2, graph = TRUE,
            main = "Mantel test: a random example with 6 X 6 matrices",
            xlab = "z-statistic", ylab = "Density",
            sub = "The vertical line shows the observed z-statistic")

# example from http://www.ats.ucla.edu/stat/r/faq/mantel_test.htm

library(ade4)
ozone<-read.table("http://www.ats.ucla.edu/stat/r/faq/ozone.csv", sep=",", header=T)

station.dists <- dist(cbind(ozone$Lon, ozone$Lat))
dim(station.dists)
ozone.dists <- dist(ozone$Av8top)

as.matrix(station.dists)[1:5, 1:5]

as.matrix(ozone.dists)[1:5, 1:5]

# 2个结果不一样
mantel.rtest(station.dists, ozone.dists, nrepet = 10000)
mantel.test(as.matrix(station.dists), as.matrix(ozone.dists),nperm = 10000)

ref:Mantel Test_王强_新浪博客 (sina.com.cn)