DALS021-MDS与PCA

作者: backup备份 | 来源:发表于2019-10-16 10:49 被阅读0次

    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最大,

    YN \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

        本文链接:https://www.haomeiwen.com/subject/cwuymctx.html