DALS021-MDS与PCA

2019-10-16 本文已影响0人 backup备份

title: DALS021-MDS与PCA
date: 2019-08-21 12:0:00
type: "tags"
tags:

MDS
PCA
categories:
Data Analysis for the life sciences

前言

这一部分是《Data Analysis for the life sciences》的第8章统计模型的第2小节，这一部分的主要内容涉及MDS和PCA，相应的Rmarkdown文档可以参考作者的Github。

MDS

MDS的全称为multi-dimensional scaling，即多维数据缩放。在这一部分中，我们会使用基因表达的数据来作为案例讲解一下。为了简化说明，我们仅考虑3个组织：

library(rafalib)
library(tissuesGeneExpression)
data(tissuesGeneExpression)
colind <- tissue%in%c("kidney","colon","liver")
mat <- e[,colind]
group <- factor(tissue[colind])
dim(mat)

结果如下所示：

> dim(mat)
[1] 22215    99

现在我们要研究一下这个数据集，我们想知道，储存在mat列中的基因表达谱的数据在不同的组织间的相似性如何。由于数据很大，无法直接画出相应的多维点图。我们通常只能绘制出二维图形，如果我们要绘制出每两个样本之间的基因表达情况不现实。而MDS图形就是为了解决这个问题而提出来的。

MDS背后的数学原理

前面我们已经知道了SVD和矩阵代数，那么我们理解MDS就相对清楚了。为了说明MDS，我们先来看一下SVD分解，如下所示：
$\mathbf{Y} = \mathbf{UDV}^\top$
我们假设 $\mathbf{U^\top Y=DV^\top}$ 的前两列的平方和剩余列的平方和。因此它们可以写为 $d_1+ d_2 \gg d_3 + \dots + d_n$ 其中 $d_i$ 是 $\mathbf{D}$ 是第i列（原文是i-th entry）。当出现这种情况时，我们就会得到如下公式：

$\mathbf{Y}\approx [\mathbf{U}_1 \mathbf{U}_2] \begin{pmatrix} d_{1}&0\\ 0&d_{2}\\ \end{pmatrix} [\mathbf{V}_1 \mathbf{V}_2]^\top$

这就表明，第 $i$ 列近似等于：
$\mathbf{Y}_i \approx [\mathbf{U}_1 \mathbf{U}_2] \begin{pmatrix} d_{1}&0\\ 0&d_{2}\\ \end{pmatrix} \begin{pmatrix} v_{i,1}\\ v_{i,2}\\ \end{pmatrix} = [\mathbf{U}_1 \mathbf{U}_2] \begin{pmatrix} d_{1} v_{i,1}\\ d_{2} v_{i,2} \end{pmatrix}$
如果我们们定义下面的二维向量：

$\mathbf{Z}_i=\begin{pmatrix} d_{1} v_{i,1}\\ d_{2} v_{i,2} \end{pmatrix}$

那么：
$\begin{align*} (\mathbf{Y}_i - \mathbf{Y}_j)^\top(\mathbf{Y}_i - \mathbf{Y}_j) &\approx \left\{ [\mathbf{U}_1 \mathbf{U}_2] (\mathbf{Z}_i-\mathbf{Z}_j) \right\}^\top \left\{[\mathbf{U}_1 \mathbf{U}_2] (\mathbf{Z}_i-\mathbf{Z}_j)\right\}\\ &= (\mathbf{Z}_i-\mathbf{Z}_j)^\top [\mathbf{U}_1 \mathbf{U}_2]^\top [\mathbf{U}_1 \mathbf{U}_2] (\mathbf{Z}_i-\mathbf{Z}_j) \\ &=(\mathbf{Z}_i-\mathbf{Z}_j)^\top(\mathbf{Z}_i-\mathbf{Z}_j)\\ &=(Z_{i,1}-Z_{j,1})^2 + (Z_{i,2}-Z_{j,2})^2 \end{align*}$
上面的这个推导告诉我们，在样本 $i$ 和样本 $j$ 之最的距离近拟等于下面二维数据点的距离：

$(\mathbf{Y}_i - \mathbf{Y}_j)^\top(\mathbf{Y}_i - \mathbf{Y}_j) \approx (Z_{i,1}-Z_{j,1})^2 + (Z_{i,2}-Z_{j,2})^2$

因为 $Z$ 是一个二维向量，因此我们可以通过绘制 $\mathbf{Z_{1}}$ 和 $\mathbf{Z_{2}}$ 来发展示这两个样本的距离。现在我们绘制出它们的距离：

s <- svd(mat-rowMeans(mat))
PC1 <- s$d[1]*s$v[,1]
PC2 <- s$d[2]*s$v[,2]
mypar(1,1)
plot(PC1,PC2,pch=21,bg=as.numeric(group))
legend("bottomright",levels(group),col=seq(along=levels(group)),pch=15,cex=1.5)

image

从图片上我们可以看出，数据点按照相应的组织区分开来了。上面的这种分开的精确近似取决于前两个主成分解释变异的程度。像上面那样所示，我们可以绘制出每个主成分可以解释的变异程度：

plot(s$d^2/sum(s$d^2))

image

虽然前两个主成分解释了超过50%的变异，不过前面的图形还是没有展示出大量的信息。但是这种图已经足够用于进行可视化大量的数据了。此外，我们还可以注意到，我们能够绘制其它的主成分来研究这些数据点，例如我们绘制第3个和第4个主成分：

PC3 <- s$d[3]*s$v[,3]
PC4 <- s$d[4]*s$v[,4]
mypar(1,1)
plot(PC3,PC4,pch=21,bg=as.numeric(group))
legend("bottomright",levels(group),col=seq(along=levels(group)),pch=15,cex=1.5)

image

从上面图形中我们可以看到，第4个主成分能够将肾脏组织的样本强烈分开。在后面的部分中，我们会讲到批次效应(batch effects)会解释这种情况。

`cmdscale()`函数

我们在上面使用了svd()函数来进行计算，不过R中有一个专门的函数用于计算MDS，生成MDS图。这个函数就是cmdscale()函数，这个函数将距离对象作为参数，然后使用主成分分析来对这些距离进行近似计算。这个函数比使用svd()函数更高效（因为不可能实现完全的svd()函数计算，那样比较花时间）。此函数默认返回二维的数据，不过我们通过设定参数k（默认情况下，k=2）可以改变结果中的维度：

d <- dist(t(mat))
mds <- cmdscale(d)
mypar()
plot(mds[,1],mds[,2],bg=as.numeric(group),pch=21,
xlab="First dimension",ylab="Second dimension")
legend("bottomleft",levels(group),col=seq(along=levels(group)),pch=15)

image

再看另外一个：

mypar(1,2)
for(i in 1:2){
plot(mds[,i],s$d[i]*s$v[,i],main=paste("PC",i))
b = ifelse( cor(mds[,i],s$v[,i]) > 0, 1, -1)
abline(0,b) ##b is 1 or -1 depending on the arbitrary sign "flip"
}

image

任意符号

SVD并非是唯一的，只要我们用-1乘以 $\mathbf{U}$ 的样本列，我们就能使用-1乘以 $\mathbf{V}$ 的任意列，通过下面的转换我们就能看出来（这一段不懂）：
$\mathbf{-1UD(-1)V}^\top = \mathbf{UDV}^\top$

扣除平均值

在所有的计算中，当我们计算SVD时，都会扣除行(row)的均值。如果我们要试图计算两列之间的近似距离，那么在 $\mathbf{Y}_{i}$ 和 $\mathbf{Y}_{j}$ 之间的距离就与 $\mathbf{Y}_i - \mathbf{\mu}$ 和 $\mathbf{Y}_j - \mathbf{\mu}$ 之间的距离相同，因为当我们过计算时，中间的 $\mu$ 就会被消去：
$\left\{ ( \mathbf{Y}_i- \mathbf{\mu} ) - ( \mathbf{Y}_j - \mathbf{\mu} ) \right\}^\top \left\{ (\mathbf{Y}_i- \mathbf{\mu}) - (\mathbf{Y}_j - \mathbf{\mu} ) \right\} = \left\{ \mathbf{Y}_i- \mathbf{Y}_j \right\}^\top \left\{ \mathbf{Y}_i - \mathbf{Y}_j \right\}$
因为扣除行均值可以降低总的变异，它可以使得SVD的结果近更为逼近。

练习

P357

PCA

PCA的相关资料可以参考作者的Github。

前面我们已经提到了PCA，这里继续深入一步，讲一下PCA背后的数学原理。

案例：双胞胎身高

我们先使用模拟数据的案例展示一个旋转，这个旋转与PCA有着很大的有关系：

library(rafalib)
library(MASS)
n <- 100
set.seed(1)
Y=t(mvrnorm(n,c(0,0), matrix(c(1,0.95,0.95,1),2,2)))
mypar()
thelim <- c(-3,3)
plot(Y[1,], Y[2,], xlab="Twin 1 (standardized height)", 
     ylab="Twin 2 (standardized height)", xlim=thelim, ylim=thelim)
points(Y[1,1:2], Y[2,1:2], col=2, pch=16)

image

这里我们专门来解释一下什么是什么成分（principla components）。

我们使用 $\mathbf{Y}$ 这个 $2 \times N$ 矩阵来表示我们的数据。这个类似于我们检测了两组基因的信息，每列表示1个样本。现在我们的任何就是，找到一个 $2 \times 1$ 向量 $\mathbf{u}_1$ ，使其满足 $\mathbf{u}_1^\top \mathbf{v}_1 = 1$ ，它能使 $(\mathbf{u}_1^\top\mathbf{Y})^\top (\mathbf{u}_1^\top\mathbf{Y})$ 最大。这个过程可以被视为每个样本，或 $\mathbf{Y}$ 向子空间 $\mathbf{u}_1$ 的投影。因此，我们需要将坐标系进行置换，使新的坐标系能够显示出最大变异。

我先试一下 $\mathbf{u}=(1,0)^\top$ 。这个投影公仅能够给出双胞胎1的身高（橘黄色）。图片标题中显示的是平方和。

mypar(1,1)
plot(t(Y), xlim=thelim, ylim=thelim,
     main=paste("Sum of squares :",round(crossprod(Y[1,]),1)))
abline(h=0)
apply(Y,2,function(y) segments(y[1],0,y[1],y[2],lty=2))
points(Y[1,],rep(0,ncol(Y)),col=2,pch=16,cex=0.75)

image

我们能否找到一个方向，使得坐标系旋转后，能够表示更高的变异？例如

$\mathbf{u} =\begin{pmatrix}1\\-1\end{pmatrix}$ 这个怎么样？它不满足 $\mathbf{u}^\top\mathbf{u}= 1$ ，因此我们可以使用另外一个向量，即
$\mathbf{u} =\begin{pmatrix}1/\sqrt{2}\\-1/\sqrt{2}\end{pmatrix}$

library(rafalib)
u <- matrix(c(1,-1)/sqrt(2),ncol=1)
w=t(u)%*%Y
mypar(1,1)
plot(t(Y),
     main=paste("Sum of squares:",round(tcrossprod(w),1)),xlim=thelim,ylim=thelim)
abline(h=0,lty=2)
abline(v=0,lty=2)
abline(0,-1,col=2)
Z = u%*%w
for(i in seq(along=w))
  segments(Z[1,i],Z[2,i],Y[1,i],Y[2,i],lty=2)
points(t(Z), col=2, pch=16, cex=0.5)

image

这个图形与双胞胎的差异有关，我们知道这个差异很少的。通常平方和我们可以确实这一点，最后我们试一下这个向量：

$\mathbf{u} =\begin{pmatrix}1/\sqrt{2}\\1/\sqrt{2}\end{pmatrix}$

u <- matrix(c(1,1)/sqrt(2),ncol=1)
w=t(u)%*%Y
mypar()
plot(t(Y), main=paste("Sum of squares:",round(tcrossprod(w),1)),
     xlim=thelim, ylim=thelim)
abline(h=0,lty=2)
abline(v=0,lty=2)
abline(0,1,col=2)
points(u%*%w, col=2, pch=16, cex=1)
Z = u%*%w
for(i in seq(along=w))
  segments(Z[1,i], Z[2,i], Y[1,i], Y[2,i], lty=2)
points(t(Z),col=2,pch=16,cex=0.5)

image

这个图形与重新缩放(re-scaled)后的平均高度有关，它有着最大的平方和。这是一个数学计算程序，它能够计算出一个 $\mathbf{v}$ ，能够使平方和最大，SVD就是这样的一个程序。

主成分

正交向量能够使平方和最大：

$(\mathbf{u}_1^\top\mathbf{Y})^\top(\mathbf{u}_1^\top\mathbf{Y})$

$\mathbf{u}_1^\top\mathbf{Y}$ 指的就是第1PC。e用于获得PC的加权(weights) $\mathbf{u}$ 指的就是因子载荷(loadings)。使用旋转这种操作，它指的就是第1PC的旋转方向。

为了获得第2PC，我们可以重复上述操作，但是残差如下：

$\mathbf{r} = \mathbf{Y} - \mathbf{u}_1^\top \mathbf{Yv}_1$

第2PC的向量含有以下性质：

$\mathbf{v}_2^\top \mathbf{v}_1=0$

它能使 $(\mathbf{rv}_2)^\top \mathbf{rv}_2$ 最大，

当 $Y$ 是 $N \times m$ 时，我们可以重复地找到第3，第4，第5，等主成分。

`prcomp`

我们已经介绍了如何使用SVD来计算PC。介理，R中有一个专门的函数可以用于找到主成分，即prcomp()，在这个案例中，数据默认中心化的，这个函数的使用如下所示：

pc <- prcomp( t(Y) )

计算出的结果与SVD相同，直到符号翻转（produces the same results as the SVD up to arbitrary sign flips，实在没理解这句话什么意思）

s <- svd( Y - rowMeans(Y) )
mypar(1,2)
for(i in 1:nrow(Y) ){
  plot(pc$x[,i], s$d[i]*s$v[,i])
}

image

因子载荷可以通过下面方式计算：

pc$rotation

计算结果如下所示：

> pc$rotation
           PC1        PC2
[1,] 0.7072304  0.7069831
[2,] 0.7069831 -0.7072304

它就相当于 (up to a sign flip？这个不懂) :

s$u

计算结果如下所示：

> s$u
           [,1]       [,2]
[1,] -0.7072304 -0.7069831
[2,] -0.7069831  0.7072304

解释的方差等价于：

pc$sdev

计算结果如下所示：

> pc$sdev
[1] 1.2542672 0.2141882

现在我们将Y转置一下，因为prcomp()函数与我们平时所用的高通量数据储存有点不太一样，平时我们的数据是列为样本，行为特征值，而prcomp()函数则是正好相反，它处理的数据列是特征值，行是样本名。

DALS021-MDS与PCA

前言

MDS

MDS背后的数学原理

`cmdscale()`函数

任意符号

扣除平均值

练习

PCA

案例：双胞胎身高

主成分

`prcomp`

猜你喜欢

热点阅读

DALS021-MDS与PCA

前言

MDS

MDS背后的数学原理

cmdscale()函数

任意符号

扣除平均值

练习

PCA

案例：双胞胎身高

主成分

prcomp

猜你喜欢

热点阅读

`cmdscale()`函数

`prcomp`