挖掘数据内部联系：相关性分析
原文链接：
相关性表示的是两个观测的数据向量之间的变化关系。一般来讲研究对象（样品或处理组）之间使用距离分析，而元素（物种或环境因子）之间进行相关性分析。两个变量之间的相关性可以用简单相关系数（例如皮尔森相关系数等）进行表示，相关系数越接近1，两个元素相关性越大，相关系数越接近0，两个元素越独立。
1. 相关系数简介
常见的相关系数及其计算方法如下所示：
①Pearson皮尔森相关系数
皮尔森相关系数也叫皮尔森积差相关系数，用来反映两个变量之间相似程度的统计量。或者说用来表示两个向量的相似度。
皮尔森相关系数计算公式如下：
分子是协方差，分母两个向量的标准差的乘积。显然是要求两个向量的标准差不为零。
当两个向量的线性关系增强时，相关系数趋于1（正相关）或者-1（负相关）。当两个变量独立时，相关系数为0。反之，不成立。比如对于Y=X2，X服从[-1,1]上的均匀分布，此时E(XY)为0，E(X)也为0，所以ρX,Y=0，但x和y明显不独立。所以“不相关”和“独立”是两回事。当Y和X服从联合正态分布时，其相互独立和不相关是等价的。
对于居中（每个数据都剪去样本均值，居中后他们的平均值就为0）的数据来说，E(X)=E(Y)=0，此时有：
即相关系数可以看作是两个随机变量的向量的夹角的cos函数。
进一步归一化X和Y向量后，||X||=||Y||=1.相关系数即为两个向量的乘积ρX,Y=X•Y
②Spearman秩相关系数
使用Pearson线性相关系数有两个局限：一是必须假设两个向量必须服从正态分布；二是取值是等距的。
对于更一般的情况有其他的一些解决方案，Spearman秩相关系数就是其中之一。Spearman秩相关系数是一种无参数（与分布无关）的检验方法，用于度量变量之间联系的强弱。在没有重复数据的情况下，如果一个变量是另一个变量的严格单调函数，则Spearman秩相关系数就是+1或者-1，称变量完全Spearman秩相关。注意这和Pearson完全相关的区别：Pearson完全相关是只有当两个变量线性关系时，Pearson相关系数为+1或者-1。
假设有以上两个样本X、Y，对X、Y样本的观察值xi，yi（i=0, 1…n）按从大到小排序，记x'i，y'i为原始xi，yi在排序后列表中的位置，x'i，y'i称为xi，yi的秩，秩次差di=x'i-y'i。不难想到，若完全正相关则di均为0，若完全负相关那么di为n+1-2i，其平方和最大，因此Spearman秩相关系数为：
此外还有Kendall秩相关系数，不再赘述。
2. 相关系数计算
计算两个数据向量或矩阵、数据框的列之间的相关性可以使用cor()函数，其使用方法如下：
cor(x, y=NULL, use="everything", method=c("pearson", "kendall", "spearman"))
其中x为向量、矩阵、数据框，若x为矩阵、数据框y可以忽略，而use为缺失值的处理方法。当x为矩阵或数据框，计算结果为元素之间的相关性矩阵。相关性矩阵对角线为1（自相关）。
此外，当具有协变量时（需要控制的干扰变量），可以使用ggm包中的pcor()函数计算偏相关系数，其使用方法如下：
pcor(u, S)
其中u为一个向量，S为变量的协方差矩阵（可以通过函数cov()计算，在这个函数里method选择相关计算方法）。若u为数值向量则第1、2个元素指定S中要计算相关性的变量下标，其余元素为协变量的下标，u也可以为字符串向量直接指定变量名字。
3. 相关系数检验
与距离不同，相关性需要进行统计检验，假如两个变量独立，那么相关系数R应该是很接近0的，那么我们认为R是服从均值为0的正态分布，那么对于实际观测值r可以构造统计量使用t检验进行分析。然而对于样本总体分布未知的时候我们计算秩相关系数，这时候最常用的方法是秩相关检验。与相关系数计算方法对应的具有相应检验方法。
在R中相关性与偏相关的检验可以通过cor.test()与pcor.test()函数分别进行，其使用方法如下所示：
cor.test(x, y,method=c("pearson","kendall","spearman"), ...)
pcor.test(r, q, n)
其中r为偏相关系数，q为协变量个数，n为样品数目。
但是这两个函数每次只能检验一个相关系数，Hmisc包中的rcorr()函数可以同时计算相关性矩阵并进行检验（具体见下一小节），同时获得相关系数矩阵与对应的p值矩阵。
4. 多重检验p值校正
假设检验是一种概率判断，因为小概率事件发生了所以我们拒绝假设。然而同时多次做这种概率判断，也会出错。例如当我们进行多重独立比较相关性时，加入有k个变量，那么需要进行k(k-1)/2个相关性分析，每个相关性均检验一次。在显著水平0.05（置信水平0.95）的情况下做出显著性判断，其正确概率为0.95，而n个独立检验均正确的概率为0.95n。若要使所有检验结果正确的概率大于0.95，则需要调整显著水平或更常用的p值校正，一个常见的方法是Bonferroni校正，其原理为在同一数据集做n个独立的假设检验，那么每一个检验的显著水平应该为只有一个检验时的1/n。例如我们只做两个变量相关检验，那么显著水平0.05，假如同时做一个数据集5个变量相关检验，因为要检验=10次，那么显著水平应为0.005，因此做Bonferroni校正后判断为显著的检验p值为原来p值的10倍。
在R中p值校正可以使用p.adjust()函数，其使用方法如下所示：
p.adjust(p, method=p.adjust.methods, n=length(p))
其中p为相关检验的结果（数值向量），n为独立检验次数，一般为length(p)，method为矫正方法，常用的方法有"bonferroni"、"holm"、"hochberg"、"hommel"、"BH"、"fdr"、"BY"、"none"。其中刚刚提到的"bonferroni"最为保守，也即校正后p值变大最多，一般不是很常用，其余方法均为各种修正方法。
校正后的p值常称为q值，使用Benjamini-Hochberg（BH）方法校正的p值也称为错误发现率（false discovery rate，FDR）。
ltm包中的rcor.test()函数在计算相关系数检验的同时还提供p值校正，其校正方法与p.adjust()函数相同，用法如下所示：
rcor.test(mat,p.adjust=FALSE, p.adjust.method="holm", ...)
其中mat为数值矩阵，p.adjust为是否需要p值校正，p.adjust.method为矫正方法。在某些很重要的多重或者多元显著性检验（例如差异基因和物种筛查）中，p值校正是必不可少的。
接下来我们以微生物群落数据为例，在R语言平台中计算物种之间以及物种与环境因子之间的Spearman相关性，并使用聚类热图进行展示，具体方法如下所示：
物种之间相关性热图
物种与环境因子相关性热图
作图代码与示例数据请关注公众号“微生态与微进化”索取。
微生态与微进化
