排序分析（Ordination Analysis） - ricket-sjtu/bi028 GitHub Wiki

1. 什么是排序分析（ordination analysis）

所谓排序分析（ordination analysis），最早是生态学（ecology）中研究群落（communities）的一大类多元分析手段，将某个地区调查的不同环境（site）以及所对应的物种组成（species），按照相似度（similarity）或距离（distance）对site在排序轴上（ordination axes）进行排序，将其表示为沿一个或多个排序轴排列的点，从而分析各个site或species与环境因子之间的关系。其目的是把多维空间压缩到低维空间（如二维），并且保证因维数降低而导致的信息量损失尽量少，实体（site或species）按其相似关系重新排列，提高其可理解性（interpretability）；同时，通过统计手段检验排序轴（ordination axes）是否能真正代表环境因子的梯度（gradient）。

如今，排序分析也被广泛应用到微生物组（microbiome）的分析中。

2. 排序分析的分类

  1. 间接排序法（indirect ordination analysis）：在特定的梯度（环境轴）上探讨物种的变化情况。
  • 主成份分析（principal component analysis, PCA）
    • stats::prcomp()
    • stats::princomp()
    • FactoMineR::PCA()
    • ade4::dudi.pca()
    • ExPosition::epPCA()
    • vegan::rda()
  • 对应分析（correspondence analysis, CA）：CA处理的是将分类数据处理成的列联表（contingency table），然后计算惯量（inertia），将惯量采用类似PCA的手段，投影到低维空间的多元统计学方法。进行CA分析的R函数有：
    • MASS::corresp()
    • FactoMineR::CA()
    • ade4::coa()
    • ca::ca()
    • ExPosition::epCA()
  • 去趋势对应分析（detrended correspondence analysis, DCA）：vegan::decorana()
  • 非度量多维缩放分析（nonmetric multidimensional scaling analysis, NMDS）
  1. 直接排序法（direct ordination analysis）：寻求潜在的（latent）或间接的环境因子来探讨物种数据的变化。
  • 冗余分析（redundancy analysis, RDA）
  • 典型对应分析（canonical correspondence analysis, CCA）
  • 去趋势典型对应分析（detrended canonical correspondence analysis, DCCA）
  • 协惯量分析（co-inertia analysis, CoIA）：分析两个数据集之间的协方差。常用的R函数包括ade4::coinertia()， made4::cia()， made4::mcoa()

其中PCA和RDA是基于线性模型（linear model）的，而CA、DCA、CCA、DCCA是基于单峰（unimodal）模型。

如何选择用单峰模型还是线性模型？

(1) 用DCA（vegan::decorana()）先对数据（site-species）进行分析；
(2) 查看结果中的“Axis lengths”的第一轴DCA1的值，根据该值判断该采用线性模型还是单峰模型：
  • 如果大于4.0，就应该选单峰模型；
  • 如果3.0-4.0之间，选线性模型或者单峰模型均可；
  • 如果小于3.0, 线性模型的结果要好于单峰模型

3. 排序分析实例

3.1 对应分析实例应用

3.2 典型对应分析实例

3.3 冗余分析实例

3.4 协惯量分析实例

这里采用的是R中的made4::cia()函数进行分析的实例：

  • made4::cia(df1, df2, cia.nf=2, cia.scan=FALSE, nsc=TRUE,...)

参数

  • df1, df2：两个类别可为“matrix”、“data.frame”、“ExpressionSet”、“marrayRaw”类型的数据集，必须有相同的行；
  • cia.nf：整数类型，协惯量分析保留的轴数量；
  • cia.scan：逻辑类型，是否绘制协惯量分析的特征值图（scree plot），以便选择cia.nf；
  • nsc：逻辑类型，是否在协惯量分析中采用非对称对应分析（Non-Symmetric Correspondence analysis, NSC）ade4::dudi.nsc()，推荐使用；如果不使用的话，则df1采用常规对应分析dudi.coa(df1)进行分析，而df2则采用行加权对应分析（Row-Weighted COrrespondence Analsyis RWCoA）ade4::rwcoa(df2, w)进行，这里的$w$为由df1对应分析返回的权重

返回结果

  • cia.result$coinertia$RV：取值在0-1之间，反映两数据集的全局相似性（global similarity）

附录：关于ordination analysis的理论知识

1. PCA

2. CA

(1) 将数据处理成列联表

CA常用于处理一组个体（individuals），多个分类变量的情况，这样我们可以构建列联表，例如对于分类变量$V_1$（类别数为$I$）和变量$V_2$（类别数为$J$），我们可以构建$I \times J$的列联表，这里每个元素$f_{ij}$分别代表对应分类$(i, j)$的频率： $$ f = \begin{pmatrix} f_{11} & f_{12} & \dots & f_{1J}\\ \vdots & \vdots & \ddots & \vdots\\ f_{I1} & f_{I2} & \dots & f_{IJ}\\ \end{pmatrix} $$

(2) 定义卡方距离

定义变量$V_1$两个类别之间的卡方距离： $$ d_{\chi^2}(i, i') = \sum_{j=1}^J \frac{1}{f_{\cdot j}}\left(\frac{f_{ij}}{f_{i\cdot}} - \frac{f_{i'j}}{f_{i'\cdot}}\right)^2; d_{\chi^2}(i, G_I) = \sum_{j=1}^J \frac{1}{f_{\cdot j}}\left(\frac{f_{ij}}{f_{i\cdot}} - f_{\cdot j}\right)^2 $$

  • 这里的距离类似欧氏距离$\left(\frac{f_{ij}}{f_{i\cdot}} - \frac{f_{i'j}}{f_{i'\cdot}}\right)^2$，但同时还考虑了一个权重，也就是每一列（column）的权重$\frac{1}{f_{\cdot j}}$；
  • $\frac{f_{ij}}{f_{i\cdot}}$相当于第$i$行在第$j$列的标准化坐标（normalized coordinate）

基于同样的原理，定义列变量$V_2$不同分类之间的距离： $$ d_{\chi^2}(j, j') = \sum_{i=1}^I \frac{1}{f_{i \cdot}}\left(\frac{f_{ij}}{f_{\cdot j}} - \frac{f_{ij'}}{f_{\cdot j'}}\right)^2 ; d_{\chi^2}(j, G_J) = \sum_{i=1}^I \frac{1}{f_{i \cdot}}\left(\frac{f_{ij}}{f_{\cdot j}} - f_{i \cdot}\right)^2 $$ 其中$G_I$和$G_J$分别为行平均（row mean profile）和列平均（column mean profile）。

如果$V_1$和$V_2$独立，则

  • $\frac{f_{ij}}{f_{i \cdot}} = f_{\cdot j}\; \forall i, j$
    • the profiles are the same as the row mean profile;
    • the profiles $N_I$ just become a single point, $G_I$, thus the cloud of the points has $\text{inertia}=0$
  • $\frac{f_{ij}}{f_{\cdot j}} = f_{i \cdot}\; \forall i, j$
    • the profiles are the same as the column mean profile;
    • the profiles $N_J$ just become a single point, $G_J$, thus the cloud of the points has $\text{inertia}=0$

(3) 计算数据与独立性假设之间的偏差

$$ \begin{array}{lcl} \text{Inertia}(N_I, G_I) & = & \sum_{i=1}^I \text{Inertia}(i, G_I)\\ & = & \sum_{i=1}^I f_{i\cdot} d_{\chi^2}(i, G_I)\\ & = & \sum_{i=1}^I \left( f_{i \cdot} \sum_{j=1}^J \frac{1}{f_{\cdot j}}\left(\frac{f_{ij}}{f_{i\cdot}} - f_{\cdot j}\right)^2 \right)\\ & = & \sum_{i=1}^I \sum_{j=1}^J \frac{(f_{ij} - f_{i\cdot} f_{\cdot j})^2}{f_{i\cdot} f_{\cdot j}}\\ & = & \frac{\chi^2}{n} = \phi^2 \end{array} $$

  • 这里$\phi^2$可以反映两个变量的关联强度（link strength）；
  • 这样惯量（$\text{Inertia}$）的计算就转化为偏离独立假设的分值（$\chi^2$）；
  • 同样的，$\text{Inertia}(N_J, G_J) = \text{Inertia}(N_I, G_I) = \phi^2$，这称为对偶（duality）；
  • 因此，我们可以说，在CA中，行和列是对称的（Symmetric）。

（4）采用因子分析对$N_I$的惯量进行分解，以最大化惯量为原则，将其投影到正交轴上

这个过程首先要确定第一个投影平面$P$，要让惯量的分量最大化，也就是需要将$N_I$中的每个点投影到该平面$P$上面的点（$H_i$）到原点（$O=G_I$）的距离$OH_i$的加权平方和最大，权重为$f_{i\cdot}$： $$ \max \; \sum_{i=1}^{I} f_{i \cdot} (OH_i)^2 $$ 确定了第一个轴$u_1$后，再同样确定第二个轴$u_2$，但必须保证$u_2 \perp u_1$。

在第$s$个正交轴上，对应的惯量值$f_{i\cdot} (OH_i^s)^2 = \lambda_s$，也就是第$s$个特征值。

（5）可视化分析（Visualization）

将$N_I$与$N_J$都分别投影到2-3个正交轴形成的平面或三维空间中，就可以进行可视化分析：

  • 该点离原点越近，表明其profile与平均值越是相似；
  • 变量1的水平$i$与变量2的水平$j$，表明$i$与$j$越是相关，反之亦然
  • 利用专业知识，分别观察不同的投影轴拆分的不同变量的不同因素

(6) 分析第$s$轴（axis）解释惯量的百分比：

$$ \frac{\text{Projected inertia of } N_I \text{ on } u_s}{\text{Total inertia of } N_I} = \frac{\sum_{i=1}^I f_{i\cdot} (OH_i^s)^2}{\sum_{i=1}^I f_{i\cdot} (OM_i)^2} = \frac{\lambda_s}{\sum_{k=1}^K \lambda_k} $$ 这里是一个例子：

       eigenvalue percentage of variance cumulative percentage of variance
dim 1 0.035401539               57.04286                          57.04286
dim 2 0.013114735               21.13190                          78.17476
dim 3 0.007301114               11.76436                          89.93912
dim 4 0.006243913               10.06088                         100.00000
  • 第一维能表示大约$57%$的独立性偏离；
  • 可以看到，偏离独立性假设的量也就是惯量可以被前面两维有效捕获（$78%$）；
  • 特征值从大到小排列，绘出对应的类似碎石图（screeplot），可用于选择需要保留的轴的数目；
  • 在CA中，特征值也就是惯量的值在0-1范围内；这跟PCA有很大的不同，在PCA中，$\lambda \ge 1$；
  • 当存在某个特征值为1，说明存在明显的结构，也就是列联表表现为分块矩阵（block），也就是说，行变量的某些类仅与列变量的某些类相关；

(7) 如何将行变量和列变量同时显示在一个图中

出于展示的要求，需将行变量和列变量同时显示在同一个坐标中，这样，需要对行、列坐标进行转换，两者符合这样的性质： $$ F_s(i) = \frac{1}{\sqrt{\lambda_s}} \sum_{j=1}^J \frac{f_{ij}}{f_{i \cdot}} G_s(j) $$ 其中

  • $F_i(s)$：变量$V_1$中，第$i$行类在第$s$轴的坐标（coordinate of row $i$ on $s$-th axis）

  • $G_j(s)$：变量$V_2$中，第$j$列类在第$s$轴的坐标（coordinate of column $j$ on $s$-th axis）

  • $\lambda_s$：第$s$个特征值（$s$-th eigenvalues, inertia associated with $s$ th axis）

  • $\frac{f_{ij}}{f_{i \cdot}}$：normalized $j$-the element of profile $i$

  • 这就相当于在第$s$轴上计算每一列$j$的重心（barycenter），权重为$\frac{f_{ij}}{f_{i\cdot}}$，然后再除以相应的$\sqrt{\lambda_s}$。

  • 这样，$\lambda_s$越小，重心就越是远离原点（origin）。

相似的， $$ G_s(j) = \frac{1}{\sqrt{\lambda_s}} \sum_{i=1}^I \frac{f_{ij}}{f_{\cdot j}} F_s(i) $$

  • 行和列的关联越是紧密，其在图中的距离就越是相近。

(8) 结果的进一步理解

$$ \frac{\text{projected inertia of } M_i \text{ on } u_s}{\text{total inertia of } M_i} = \frac{f_{i\cdot} (OH_i^s)^2}{f_{i\cdot} (OM_i)^2} = \cos^2(\overrightarrow{OM_i}, u_s) $$

  • $\overrightarrow{OM_i}$：从坐标原点$O$指向点$M_i$的向量
  • $OH_i^s$：$M_i$投影到第$s$轴的点$H_i$到坐标原点的距离
  • $\cos(\cdot, \cdot)$：两个向量夹角的余弦（cosine）

实例代码

library(FactoMineR)
data(children)
ca.children <- CA(children, row.sup=15:18, col.sup=6:8)
summary(ca.children)
chi2 <- sum(ca.children$eig[,1]) * sum(children[1:14,1:5], na.rm=TRUE) # \chi^2 = \phi * n, n=#obs
pval <- 1 - pchisq(chi2, df=(14-1)*(5-1)) # chi^2 dist with df=(nrow-1)(ncol-1) for contingency table
cat("The row variable has 14 categories, while the column variable has 5 categories; \n>and the deviation from the indpependence is chi^2=", chi2, "\n>and the corresponding p-value=", pval, " with df=52\n")
ellipseCA(ca.children, 
          ellipse="col", # show column ellipse
          col.col.ell=c(rep("blue",2),rep("transparent",3)), # set color for ellipses
          invisible=c("row.sup","col.sup") # set some variables invisible
)

另外一个自己构建的实例：

test <- data.frame(matrix(c(10,0,0,0,9,1,0,3,7), nrow=3, byrow=T), 
                 row.names=c("Sweet", "Sour", "Bitter"))
colnames(test) <- c("Perc.sweet", "Perc.sour", "Perc.bitter")
ca.test <- CA(test)

我们可以看到

  • 在第一维，可以很容易将sweet与其他sour、bitter分开，因为其eigenvalue是1；
  • 在第二维，可以将sour和bitter分开，当两者联系加强后，两者在第二维的距离变近，导致判别能力变弱，对应的eigenvalue也减小；
  • 如果两个变量具有非常强的相关性，则总的inertia会接近保留轴axes的数量n.axes；
  • 这里有排他性相关（exclusive association）的概念，也就是这里sweet只会被识别为”perceived sweet“。 这就是一个分块的例子，我们可以将sour与bitter的联系加强：
test <- data.frame(matrix(c(10,0,0,0,7,3,0,4,6), nrow=3, byrow=T), 
                 row.names=c("Sweet", "Sour", "Bitter"))
colnames(test) <- c("Perc.sweet", "Perc.sour", "Perc.bitter")
ca.test <- CA(test)

另一种极端的情况（对角矩阵）：

test <- data.frame(matrix(c(10,0,0,0,9,0,0,0,8), nrow=3, byrow=T), 
                 row.names=c("Sweet", "Sour", "Bitter"))
colnames(test) <- c("Perc.sweet", "Perc.sour", "Perc.bitter")
ca.test <- CA(test)

观察一下点图在Dim1和Dim2发生了什么变化？

3. MCA (Multiple correspondence analysis)

MCA与CA不同，其应对的数据是下面的这种数据： $$ f = \begin{pmatrix} f_{11} & f_{12} & \dots & f_{1J}\\ \vdots & \vdots & \ddots & \vdots\\ f_{I1} & f_{I2} & \dots & f_{IJ}\\ \end{pmatrix} $$ 包含了$I$个个体（individuals），$J$个分类变量，其中第$j$个分类变量具有$K_j$个类，需要将对应的数据进行转换。

(1) 转换为指示矩阵

$$ y = \begin{pmatrix} y_{11} & y_{12} & \dots & y_{1K}\\ \vdots & \vdots & \ddots & \vdots\\ y_{I1} & y_{I2} & \dots & y_{IK}\\ \end{pmatrix} $$

  • 假设对于第$j$个变量，假设有4个分类，则需要将其转换为“1 0 0 0”、“0 1 0 0”……这种分类指示向量；
  • 这样矩阵的列数就变为$K=\sum_{j=1}^J K_j$，源矩阵就变为$I \times K$的矩阵