R语言 | 相关性分析
推荐阅读
相关性是指两个变量的关联程度，定量变量之间的关系可以用相关系数来描述。相关系数的符号（±）表明关系的方向（正相关或负相关），其值的大小表示关系的强弱程度（完全不相关时为0，完全相关时绝对值为1）
本文以某微生物组数据为例，通过在R中计算微生物-微生物相关性或微生物-环境因子的相关性，简介相关性分析在R中的运行过程。
示例数据、R脚本等，已上传至百度盘（提取码ouor）：
https://pan.baidu.com/s/1vspGJlOkX8zreKjAE1fS8A
示例文件简要
示例数据文件“data.txt”的内容展示如下。
sample，样本名称；
SOC、TN、pH等，各样本中的环境因子数据，示例文件中展示了6种环境因子数据；
Richness，各样本中，alpha多样性中的richness指数；
Proteobacteria、Acidobacteria等，各样本中的主要微生物类群，示例文件中展示了10种微生物门类群数据（界门纲目科属种的门水平）。
接下来，我们在R中运行相关性分析，计算微生物-微生物之间的相关性，以及环境因子-微生物之间的相关性等。
在R中运行相关性分析示例
R作为一门强大的统计学语言，可以计算多种相关系数，例如Pearson相关系数、Spearman相关系数、Kendall相关系数、偏相关系数、多分格（Polychoric）相关系数、多系列（Polyserial）相关系数等等。
以下简介几种常见相关系数的计算命令。
首先将示例数据读入R中。为了方便后续的演示，除了混合数据（包含环境因子数据及微生物类群丰度数据）外，再分别将环境因子数据以及细菌类群丰度数据单独提取出并赋值为新数据框。
#读取数据dat <- read.table('data.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)dat_env <- dat[1:6] #只包含环境因子数据dat_phylum <- dat[8:17] #只包含细菌类群丰度数据
cov()计算协方差（Pearson、Spearman、Kendall）
在R中，可使用cov()计算协方差，得到协方差矩阵。
例如使用cov()分别计算Pearson、Spearman、Kendall相关的协方差矩阵。关于cov()的详情可使用?cov()查看帮助。
#协方差计算，cov()cov_pearson <- cov(dat, method = 'pearson')cov_spearman <- cov(dat, method = 'spearman')cov_kendall <- cov(dat, method = 'kendall')存储了两两变量之间的协方差数据结果。
cor()计算相关性（Pearson、Spearman、Kendall）
与cov()命令类似，使用cor()命令即可计算Pearson、Spearman、Kendall相关系数。关于cor()的详情可使用?cor()查看帮助。
#相关系数计算，cor()cor_pearson <- cor(dat, method = 'pearson')cor_spearman <- cor(dat, method = 'spearman')cor_kendall <- cor(dat, method = 'kendall')存储了两两变量之间的Spearman相关系数，各变量之间的相关性清晰可见。
在得到了相关系数结果后，若有需要，我们可以从中做一些筛选，并将结果输出在本地。
#例如提取“cor_pearson”中相关系数 >0.5 或 <-0.5 的结果输出为 csv 样式cor_pearson[abs(cor_pearson) <= 0.5] <- 0write.csv(cor_pearson, 'cor_pearson.csv', quote = FALSE)正常情况下，相关系数太低的数据一般不是我们想要的结果，可以过滤掉。此处选择保留Pearson相关系数绝对值大于等于0.5的值，并将绝对值低于0.5的均设为0后，输出结果如下所示。我们即可根据筛选后的相关性分析结果查看所关注的重要变量之间是否存在较强的相关性，以进行后续的统计分析工作。
注：作为示例，此处的筛选方法仅供参考，大家在实际的数据分析中可自行灵活决定过滤标准。
如上所述，默认情况下使用cor()命令处理单一数据框时，得到的结果是一个对称矩阵，数据框中所有变量之间两两计算相关。
我们还可以通过在cor()命令中输入两个包含不同变量的数据框，计算两个数据框中变量相互之间的相关系数，得到非对称矩阵。如下示例，使用cor()计算上文数据框“dat_phylum”（只包含细菌类群丰度数据）和“dat_env”（只包含环境因子数据）中数据间的相关性，即此处着重关注细菌类群丰度与环境因子之间的相关性。
#指定分组的相关性分析#此处计算“dat_phylum”和“dat_env”中数据的相关性，即只关注细菌类群丰度与环境因子之间的相关性phylum_env_spearman <- cor(dat_phylum, dat_env, method = 'spearman')查看所得到的结果“phylum_env_spearman”，如下所示。该结果不再为一个对称矩阵，只包含了我们所期望得到的微生物-环境因子间的相关系数。同样地，若有需要，可对相关系数作一定筛选后write.table()输出在本地以供后续分析等，不再多说。
pcor()计算偏相关
偏相关是指在控制一个或多个定量变量时，另外两个定量变量之间的相互关系。R包ggm中提供的命令pcor()可以计算偏相关系数。
该命令调用格式为：“pcor(u, S)”
其中，u是一个向量，向量中前两个元素为要计算相关系数的变量下标（或名称），其余元素为条件变量（即要排除影响的变量）的下标（或名称）；S为变量的协方差矩阵。可在加载ggm包后使用?pcor()查看命令详情。
以下示例使用pcor()计算上文数据框“dat_phylum”（只包含细菌类群丰度数据）中，Proteobacteria和Acidobacteria的偏相关系数。
#此处计算“dat_phylum”中，Proteobacteria 和 Acidobacteria 的偏相关系数library(ggm)
select <- c('Proteobacteria', 'Acidobacteria')delet <- names(dat_phylum)[-which(names(dat_phylum) %in% select)]pcor(c(select, delet), cov(dat_phylum, method = 'spearman'))计算结果屏幕输出如下。即在控制了其它微生物类群的影响后，Proteobacteria和Acidobacteria的相关系数约为-0.812，二者存在较强的负相关关系。
其它类型的相关
除了常见的Pearson、Spearman、Kendall等相关系数，更多类型的相关系数同样可以在R中计算得到。
例如，polycor包中的hetcor()函数可以计算一种混合的相关矩阵，其中包括数值型变量的Pearson积差相关系数、数值型变量和有序变量之间的多系列相关系数、有序变量之间的多分格相关系数以及二分变量之间的四分相关系数。多系列、多分格和四分相关系数都假设有序变量或二分变量由潜在的正态分布导出。
关于hetcor()此处不对其作详细介绍，若有需要可使用?hetcor()参阅其R文档。
以下示例为使用hetcor()计算上文数据框“dat_phylum”（只包含细菌类群丰度数据）中各细菌类群之间的多分格（Polychoric）相关系数的一个简单示例。
#hetcor()，以计算“dat_phylum”中两两细菌类群间多分格相关系数为例library(polycor)
cor_polychoric <- hetcor(dat_phylum, type = 'Polychoric')cor_polychoric$correlations计算多分格相关系数后，将计算结果存储在变量“cor_polychoric”中。其中两两细菌类群间的多分格相关系数值即可使用“cor_polychoric$correlations”查看。同样地，若有需要可将相关系数矩阵提取后write.table()输出在本地以供后续分析使用。
相关性的显著性检验
上文中简介了几种相关系数在R中的计算方法。
严格来说，在计算好相关系数后，还应当进行统计显著性检验，以判断变量之间的相关性是否显著，排除假阳性结果。对于相关系数的显著性检验同样存在多种方法（命令）可供选择。
cor.test()对相关系数执行检验
cor.test()可对单个的Pearson、Spearman、Kendall相关系数进行检验，使用?cor.test()查看命令详情。
如下示例，使用cor.test()判断上文计算得到的相关系数矩阵“cor_pearson”中Proteobacteria和Acidobacteria的相关系数是否显著。
#此处使用“dat”中 Proteobacteria 和 Acidobacteria 列的数据#意在判断“cor_pearson”中 Proteobacteria 和 Acidobacteria 的相关系数是否显著cor.test(dat[,'Proteobacteria'], dat[,'Acidobacteria'], method = 'pearson')在结果中我们观察到p值为6.686e-07远低于显著性水平（默认为0.05，零假设变量间不相关，p值显著即可拒绝零假设），因此可判断Proteobacteria和Acidobacteria的Pearson相关系数显著。
corr.test()计算相关性并执行检验
如上命令cor.test()一次只能检验一组数据间的相关关系，数据较多时使用起来并不方便（即便是以循环的方式处理）。
R包psych中的命令corr.test()提供了在计算变量间相关系数的同时执行显著性检验的方法，可一次处理完所有的组。该命令的用法类似于cor()，更多详情可使用?corr.test()查看。
如下示例，分别使用corr.test()处理单个数据框，计算上文数据框“dat_phylum”（只包含细菌类群丰度数据）中各细菌类群之间的Spearman相关性，得到包含两两细菌类群之间Spearman相关系数的对称矩阵；以及处理两个独立的数据框，计算“dat_phylum”（只包含细菌类群丰度数据）和“dat_env”（只包含环境因子数据）中数据间的Spearman相关性，即计算细菌类群-环境因子之间的Spearman相关系数。
#相关性的显著性检验library(psych)
#计算phylum_corr <- corr.test(dat_phylum, method = 'spearman')phylum_env_corr <- corr.test(dat_phylum, dat_env, method = 'spearman')
#以“phylum_corr”为例，查看数据phylum_corr$r #相关系数值phylum_corr$p #显著性 p 值计算结果中既包含了相关系数，又包含了显著性检验的p值。默认情况下，对于某两个变量间的相关系数，若其检验p值低于0.05则表明相关性显著（零假设变量间不相关，p值显著即可拒绝零假设）。我们可根据结果中显著性检验的p值对相关系数进行筛选，仅保留显著的相关系数以排除假阳性。
注：若数据量较大不便观测的话，可将计算结果的相关系数以及显著性p值分别赋值为数据框的类型后write.table()输出在本地查看。
以下为一个简要的筛选，以计算结果“phylum_corr”为例，根据显著性检验的p值（仅保留p<0.05的相关系数）以及相关性r值（即根据相关性强度筛选，仅保留绝对值大于等于0.5的相关系数）对相关系数做取舍。满足条件的保留，不满足条件的将转化为0，并将最终的相关系数矩阵输出在本地。
#筛选，例如在“phylum_corr”中，根据 p(<0.05) 值和 r(>=0.3 or <=-0.3) 值做保留phylum_corr$p[phylum_corr$p >= 0.05] <- -1phylum_corr$p[phylum_corr$p < 0.05 & phylum_corr$p >= 0] <- 1phylum_corr$p[phylum_corr$p == -1] <- 0phylum_corr$r[abs(phylum_corr$r) < 0.3] <- 0phylum_corr_final <- phylum_corr$r * phylum_corr$pwrite.csv(phylum_corr_final, 'phylum_corr_final.csv', quote = FALSE)最终所得相关系数矩阵“phylum_corr_final”如下所示。与一开始得到的相关系数矩阵相比，不显著的相关性以及较弱的相关性均被剃除。
注：作为示例，此处的筛选方法仅供参考，大家在实际的数据分析中可自行灵活决定过滤标准。
相关性的可视化示例（corrplot()）
若有需要，我们可选将变量间的相关性分析结果作图展示出来，以直观地观测。
在相关图的绘制方面，R中同样存在很多方法（命令）可供使用，此处不再作更多介绍了，大家可自行学习。例如，我在“R语言中相关图的绘制方法简介”中对几个常用的相关图作图R命令做了简介。
例如，我本人习惯使用corrplot包中的corrplot()命令对相关性数据进行可视化展示。如下简要展示两例corrplot()命令绘制相关图的示例，作图数据分别对应上述的结果“phylum_corr_final”以及“phylum_env_corr”。
##相关图绘制示例，corrplot 包library(corrplot)
#第一个示例，根据筛选后的数据“phylum_corr_final”作图corrplot(phylum_corr_final, method = 'number', number.cex = 0.8, diag = FALSE, tl.cex = 0.8)corrplot(phylum_corr_final, add = TRUE, type = 'upper', method = 'pie', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
#第二个示例，在作图时，通过指定过滤参数（p<0.05），作图时自动屏蔽不显著的值corrplot(phylum_env_corr$r, p.mat = phylum_env_corr$p, method = 'number', insig = 'blank', sig.level = 0.05, addCoef.col = 'black', number.cex = 0.8, tl.cex = 0.8)第一张示例图，作图数据为上文中的“phylum_corr_final”，即根据显著性检验的p值（仅保留p<0.05的相关系数）以及相关性r值（即根据相关性强度筛选，仅保留绝对值大于等于0.5的相关系数）对相关系数做取舍后的相关矩阵。根据数据内容直接读取作图。
第二张示例图，作图数据为上文中的“phylum_env_corr”。该数据中预先并未依据p值或r值作筛选，但由于corrplot()命令中自带根据p值的过滤功能，我们在作图时启动了该功能，因此在作图结果中屏蔽了p值不显著的相关系数数据。
我就知道你在看👇