挖掘数据内部联系：相关性分析

相关系数简介

常见的相关系数及其计算方法如下所示：

①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秩相关系数，不再赘述。

相关系数计算

计算两个数据向量或矩阵、数据框的列之间的相关性可以使用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也可以为字符串向量直接指定变量名字。

相关系数检验

与距离不同，相关性需要进行统计检验，假如两个变量独立，那么相关系数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值矩阵。

多重检验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相关性，并使用聚类热图进行展示，具体方法如下所示：

#读取物种和环境因子数据（行名字均是样品名，否则需要转置）
community=read.table(file="taxonomy.txt", header=T, check.names=FALSE)
rownames(community)=community[,1]
com=community[,-1]
environment=read.table(file="environment.txt", header=T, check.names=FALSE)
rownames(environment)=environment[,1]
env=environment[,-1]
#整合数据并计算相关性
data=data.frame(com, env)
data=as.matrix(data)
library(Hmisc)
corr=rcorr(data, type="spearman")
rcorr=corr$r
pcorr=corr$P
#一、抽取物种之间相关性及其p值作图
n=length(colnames(com))
spcor=as.matrix(rcorr[1:n, 1:n])
spcop=as.matrix(pcorr[1:n, 1:n])
library(gplots)
heatmap.2(spcor, col=bluered(75), srtCol=45, trace="none", key.title=NA, dendrogram="both", cexRow=1.2, cexCol=1.2, keysize=0.9, key.ylab=NA, margins=c(8, 15), key.xlab="Correlation", offsetRow=-0.1, offsetCol=0.1, cellnote=round(spcop,2), notecol='black', notecex=1, density.info="none")
#其中cellnote在热图色块中显示p值，round(x,2)表示保留x小数点后两位
#二、抽取物种与环境因子的相关性及其p值进行作图
m=length(colnames(env))
ecocor=as.matrix(rcorr[1:n, (n+1):(n+m)])
ecocop=as.matrix(pcorr[1:n, (n+1):(n+m)])
#接下来将p值用显著性符号表示
sigcor=ecocop
sigcor[which(sigcor<0.001)]="***"
sigcor[which(sigcor>=0.001&sigcor<0.01)]="**"
sigcor[which(sigcor>=0.01&sigcor<0.05)]="*"
sigcor[which(sigcor>=0.05&sigcor<0.1)]="."
sigcor[which(sigcor>=0.1)]=" "
#其中which也可以省略：sigcor[sigcor<0.001]="***"
par(mar=c(2,2,2,2))
heatmap.2(ecocor, col=bluered(75), srtCol=35, trace="none", key.title=NA, dendrogram="both", cexRow=1.2, cexCol=1.2, keysize=0.9, key.ylab=NA, margins=c(8, 15), key.xlab="Correlation", offsetRow=-0.1, offsetCol=-0.1, cellnote=sigcor, notecol='black', notecex=2, density.info="none")
#添加尾注
text(0.35, 0.035,labels="Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1")

其中物种之间相关性热图如下所示：

物种与环境因子之间的相关性热图如下所示：

示例数据链接：

https://pan.baidu.com/s/1YWwfAs6i8xV8YJzLmqkYYQ

提取码：xcfx