R包vegan的典范相关分析(CCorA)
然后计算降维后数据集之间的典范相关系数。
即可得到典范相关系数。如果上述降维为多个正交轴(最大min[p,q]),则得到各轴的典范相关系数。
CCorA(canonical correlation analysis)出现时间较早,早期也称为CCA。但考虑到后来出现了典范对应分析(canonical correspondence analysis,CCA),为了区分二者,尽量避免使用CCA作为CCorA的简称。
CCorA属于对称分析(symmetrical analysis),两组数据集相互之间无响应变量和解释变量之分,二者地位等同。例如基因共表达、物种间互作、群落演替过程中物种-环境长期相互作用等。因此,CCorA是一种用于描述数据的探索性分析方法,而非确立因果关系的模型,不涉及假设检验过程。但如果有必要,仍可以通过置换检验确立相关性的显著性。
如果分析的两个数据集之间存在因果关系,例如常见的“处理-响应”类型,或者如物种(代表响应变量)与环境(代表解释变量)之间的关系,推荐使用RDA、CCA等模型方法。
CCorA还可用于确定一组典范变量,即每个组内变量的正交线性组合,可以最好地解释组内和组间的变异性。
CCorA要求两组数据为符合多元正态分布的定量数据。且数据集中的变量具有不同的量纲时,需要标准化处理。
R包vegan的CCorA
本篇简介R包vegan的CCorA方法。
library(vegan)
#查看文档,最下方的示例
?CCorA
数据集
vegan包的内置数据集“mite”,记录了70个土壤观测地点的35种甲螨的个体数量,行为观测地点,列为甲螨物种。
library(vegan)
#甲螨数据集
#详情 https://www.rdocumentation.org/packages/vegan/versions/2.4-2/topics/mite
data(mite)
head(mite[ ,1:10])
接下来使用该甲螨数据,展示vegan中的CCorA过程。
按照文档中的示例,应该是期望确立两组亚群甲螨物种之间是否存在某种关联(正相关可表示趋于共存,负相关可表示趋于竞争)。且两组数据均为物种变量,地位等同。
#对土壤地点设置分组(分组参考自 Legendre 2005, Fig. 4 的 PCA)group.1 <- c(1, 2, 4:8, 10:15, 17, 19:22, 24, 26:30)
group.2 <- c(3, 9, 16, 18, 23, 25, 31:35)
#按分组将物种多度数据拆分为两组数据集
#物种多度数据分析,推荐执行 Hellinger 转化
mite.hel.1 <- decostand(mite[,group.1], 'hel')
mite.hel.2 <- decostand(mite[,group.2], 'hel')
rownames(mite.hel.1) = paste('S', 1:nrow(mite), sep = '')
rownames(mite.hel.2) = paste('S', 1:nrow(mite), sep = '')
多元正态性评估
CCorA前,需确定两组定量数据是否符合多元正态分布,可通过前述MANOVA中提到的Q-Q图作评估。
#Q-Q 图评估多元正态性
#若数据服从多元正态性,则点将落在直线附近
par(mfrow = c(1, 2))
qqplot(qchisq(ppoints(nrow(mite.hel.1)), df = ncol(mite.hel.1)), mahalanobis(mite.hel.1, colMeans(mite.hel.1), cov(mite.hel.1)))
abline(a = 0, b = 1)
qqplot(qchisq(ppoints(nrow(mite.hel.2)), df = ncol(mite.hel.2)), mahalanobis(mite.hel.2, colMeans(mite.hel.2), cov(mite.hel.2)))
abline(a = 0, b = 1)
“mite.hel.2”符合,但“mite.hel.1”偏离程度比较明显。
理论上该“mite.hel.1”数据集不能继续用于CCorA,但由于本文仅为示例演示,所以请允许我们忽略多元正态性的问题,继续演示。
实际的情况中,不妨试下其它的转化数据方法可不可行,否则不可使用。
CCorA执行及结果概述
vegan中,CCorA()用于执行CCorA。
#CCorA,详情 ?CCorA
#Y 和 X 地位等同
#这里两组均为物种丰度,量纲一致,且 Hellinger 转化降低了高丰度物种的权重,故无需对两组数据集标准化
#若量纲不同(比如不同的环境变量),则需要将对应的数据集标准化
#999 次置换确定相关性的显著性
out <- CCorA(Y = mite.hel.1, X = mite.hel.2, stand.Y = FALSE, stand.X = FALSE, permutations = 999)
out
Rank显示了Y和X矩阵中的变量个数。
Pillai's trace是典范相关系数的和。
Significance of Pillai's trace为置换检验的结果,其中from F-distribution代表了随机置换过程中所得的期望Pillai's trace,based on permutations代表了p值,这里可知相关性是显著的。
Canonical Correlations代表了各个轴的典范相关系数,它们的加和即等于Pillai's trace。通过前几个有代表性的轴所示的相关性,两组数据集整体相关程度很高。
RDA R squares,Y|X和X|Y的RDA的双变量冗余系数(R2),adj. RDA R squares是校正后的R2。它们严格来说不属于CCorA的计算部分,但可以评估一组变量的总变差能够由另一组变量在一个或多个正交轴上的解释程度。因为即使R2很小时,也可能出现较大的典范相关系数,此时可根据R2评估典范相关系数的合理性。
对于主要信息的提取。
#根据 summary() 的提示,提取主要信息summary(out)
#例如
#Pillai's trace
out$Pillai
#各个轴的典范相关系数
out$Eigenvalues
#校正后的 RDA R2
out$RDA.adj.Rsq
#置换检验 p 值
out$p.perm
关于理解变量间相关性
上面给出了两个数据集整体的相关性。
对于各变量间的相关性,例如我们想知道哪些物种和哪些物种之间存在较大的相关,可通过排序图观测。
vegan中提供了一些可视化方案,例如biplot()。biplot()一次会输出两个图,分别代表矩阵Y和X。由于CCorA是对称分析,所以两个排序图互有关联,典范轴显示了两个矩阵的共同变化趋势。每个排序图都是两组变量线性组合的结果,所以单个变量解读比较困难,需两个图结合起来解读。
#作图观测,如 biplot(),详情 ?CCorA
#“objects”生成两个对象图,第一个表示 Y 矩阵中的对象,第二个表示 X 矩阵中的对象;以展示第 1、2 轴为例
biplot(out, plot.type = 'objects', plot.axes = c(1, 2))
biplot(out, plot.type = 'variables', plot.axes = c(1, 2), cex = c(0.7, 0.6))
biplot(out, plot.type = 'ov', plot.axes = c(1, 2), cex = c(0.7, 0.6))
biplot(out, plot.type = 'biplots', plot.axes = c(1, 2), cex = c(0.7, 0.6))
好了我们结合上述这两个双序图作个解读,加深对这种变量间相关性的理解。CCorA中,变量间的相关性是通过在两个数据集中,变量在对象中的相似分布特征,作为呈现的。因此,排序图中的相关性可以这样解读。
参考资料