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分解,如下所示:
我们假设 的前两列的平方和剩余列的平方和。因此它们可以写为 其中 是 是第i列(原文是i-th entry)。当出现这种情况时,我们就会得到如下公式:
这就表明,第列近似等于:
如果我们们定义下面的二维向量:
那么:
上面的这个推导告诉我们,在样本和样本之最的距离近拟等于下面二维数据点的距离:
因为是一个二维向量,因此我们可以通过绘制和来发展示这两个样本的距离。现在我们绘制出它们的距离:
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
乘以的样本列,我们就能使用-1
乘以的任意列,通过下面的转换我们就能看出来(这一段不懂):
扣除平均值
在所有的计算中,当我们计算SVD时,都会扣除行(row)的均值。如果我们要试图计算两列之间的近似距离,那么在和之间的距离就与和之间的距离相同,因为当我们过计算时,中间的就会被消去:
因为扣除行均值可以降低总的变异,它可以使得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)。
我们使用 这个 矩阵来表示我们的数据。这个类似于我们检测了两组基因的信息,每列表示1个样本。现在我们的任何就是,找到一个 向量 ,使其满足 ,它能使 最大。这个过程可以被视为每个样本,或向子空间 的投影。因此,我们需要将坐标系进行置换,使新的坐标系能够显示出最大变异。
我先试一下 。这个投影公仅能够给出双胞胎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
我们能否找到一个方向,使得坐标系旋转后,能够表示更高的变异?例如
这个怎么样?它不满足 ,因此我们可以使用另外一个向量,即
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
这个图形与双胞胎的差异有关,我们知道这个差异很少的。通常平方和我们可以确实这一点,最后我们试一下这个向量:
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)后的平均高度有关,它有着最大的平方和。这是一个数学计算程序,它能够计算出一个 ,能够使平方和最大,SVD就是这样的一个程序。
主成分
正交向量能够使平方和最大:
指的就是第1PC。e用于获得PC的加权(weights) 指的就是因子载荷(loadings)。使用旋转这种操作,它指的就是第1PC的旋转方向。
为了获得第2PC,我们可以重复上述操作,但是残差如下:
第2PC的向量含有以下性质:
它能使 最大,
当 是 时,我们可以重复地找到第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()
函数则是正好相反,它处理的数据列是特征值,行是样本名。
网友评论