读书笔记（8）降维

Distance and Dimension Reduction

这是最近书上学习降维部分时的读书笔记

每个实数矩阵都有一个奇异值分解，但不一定都有特征分解。例如，非方阵的矩阵没有特征分解，这是我们只能使用奇异值分解。

高维数据的distance

下面是结果是相同的

x <- e[,1]
y <- e[,2]
sqrt(sum((x-y)^2))
sqrt( crossprod(x-y) )

d <- dist(t(e)) #“dist”
#因为dist计算行与行间距离，这里要计算样本间距离，所以转置
as.matrix(d)[1,2] #要使用行和列访问，需要强制转换为矩阵

降维动机

降维是想要将数据集缩小为较少的维度，同时大致保留重要的属性，例如样本之间的距离。如果我们能够缩小到二维，那么我们可以轻松地绘制图。

NOTE:Rotations（旋转）之后（z=Ay），每个点的standard deviations会被缩放，比如，假设y的每一列为具有标准偏差σ的独立随机变量$$sd[z_1]=sd[(y_1+y_2)/2]=\frac{1}{\sqrt{2}}\sigma$$$$sd[z_2]=sd[y_1-y_2]=\sqrt{2}\sigma$$需要进行缩放，只要A为正交矩阵$$A^{T}A=E$$因为$$A^{-1}A=E$$，所以$$A^T=A^{-1}$$它可以保证SD保留属性。

SVD（Singular Value Decomposition）

$$U^{T}Y=D{V}^T$$U为m*p正交矩阵，Y为m*n矩阵，D为p*p对角矩阵，V为n*p正交矩阵，p=min(m,n)。
t(U)为data Y的rotation，且可以使得$$U^{T}Y=D{V}^T$$的列的variability（sum of squares）是降低的。也可以写为YV=UD，其列的variability（sum of squares）也是降低的。

s <- svd(Y)
U <- s$u
V <- s$v
D <- diag(s$d) 
plot(s$d)

the sum of squares of UD 的后面几个近似于0，用前面几个就可以近似。

Yhat <- U %*% D %*% t(V)
resid <- Y - Yhat
max(abs(resid))
## [1] 3.508305e-14
k <- ncol(U)-4
Yhat <- U[,1:k] %*% D[1:k,1:k] %*% t(V[,1:k])
resid <- Y - Yhat 
max(abs(resid))
## [1] 3.508305e-14

每一列解释的variability

plot(s$d^2/sum(s$d^2)*100,ylab="Percent variability explained")
plot(cumsum(s$d^2)/sum(s$d^2)*100,ylab="Percent variability explained",ylim=c(0,100),。type="l")
降维近似未能解释的变异，下面结果一样
```javascript
var(as.vector(resid))/var(as.vector(Y))
1-sum(s$d[1:k]^2)/sum(s$d^2)

D也表明每个PC在解释可变性方面的贡献。如果y的两列相关性很高，二者的平均值更有用。

坐标旋转coordinate rotation

坐标（coordinates）是新basis定义的空间的投影。$$A^{-1}Y=Z$$Z和Y 携带相同的信息，但坐标系不同。（在线性代数中，矩阵与线性变换存在一一对应的关系，将坐标与一个矩阵相乘，相当于进行投影或旋转）

原始数据为这个样子，想要对它进行降维，可以旋转，使得变异尽量由第一维呈现。

旋转后

Multidimensional scaling

$$\mathbf{Y}={\mathbf{U}\mathbf{D}\mathbf{V}}^{\top }$$考虑用前两个主成分来近似：$$\mathbf{Y}\approx [{\mathbf{U}}_1{\mathbf{U}}_2]\left(\begin{array}{cc}{d}_1& 0\\ 0& d_2\end{array}\right)[{\mathbf{V}}_1{\mathbf{V}}_2{]}^{\top }$$Y中的第 i列$${\mathbf{Y}}_i\approx [{\mathbf{U}}_1{\mathbf{U}}_2]\left(\begin{array}{cc}{d}_1& 0\\ 0& d_2\end{array}\right)\left(\begin{array}{c}{v}_{i,1}\\ v_{i,2}\end{array}\right)=[{\mathbf{U}}_1{\mathbf{U}}_2]\left(\begin{array}{c}{d}_1{v}_{i,1}\\ d_2{v}_{i,2}\end{array}\right)$$二维向量为$${\mathbf{Z}}_i=\left(\begin{array}{c}{d}_1{v}_{i,1}\\ d_2{v}_{i,2}\end{array}\right)$$

表明两个样本间的距离可以用两个样本点之间的距离来近似。$$({\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$$

专门为MDS图制作的特殊函数。 它以距离对象作为参数，然后使用主成分分析来提供可以使用k个维度获得的最佳距离近似值。 此功能效率更高，因为不必执行完整的SVD，这很耗时。 默认情况下，它返回二维，但是我们可以通过默认为2的参数k进行更改。

d <- dist(t(mat))
mds <- cmdscale(d)
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)

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

SVD不是唯一的，因为我们可以将V的任何列乘以-1，只要我们将U的样本列乘以-1。$${-1\mathbf{U}\mathbf{D}(-1)\mathbf{V}}^{\top }={\mathbf{U}\mathbf{D}\mathbf{V}}^{\top }$$
计算两个样本间的距离，每行减去均值，再进行SVD不会有影响，因为会相互抵消。但删除行平均值会减少总variation，可以使SVD近似得更好。 { ( Y i − μ ) − ( Y j − μ ) } ⊤ { ( Y i − μ ) − ( Y j − μ ) } = { Y i − Y j } ⊤ { Y i − Y j } \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\} {(Yi​−μ)−(Yj​−μ)}⊤{(Yi​−μ)−(Yj​−μ)}={Yi​−Yj​}⊤{Yi​−Yj​}

Principal Component Analysis

计算主成分的方法：
$${\mathbf{u}}_1^{\top }\mathbf{Y}$$为第一主成分，u为loadings或称为rotations，它表示的是第一主成分的方向，第一主成分就是新的坐标。
首先，寻找使平方和最大化的正交向量：
$$({\mathbf{u}}_1^{\top }\mathbf{Y}{)}^{\top }({\mathbf{u}}_1^{\top }\mathbf{Y})$$再使用残差（residuals）$$\mathbf{r}=\mathbf{Y}-{\mathbf{u}}_1^{\top }{\mathbf{Y}\mathbf{v}}_1$$再最大化$$({\mathbf{r}\mathbf{v}}_2{)}^{\top }{\mathbf{r}\mathbf{v}}_2$$其中，v2为与v1正交的向量，且为单位向量。

下面两个结果相同（忽略符号）

pc <- prcomp( t(Y) )#行为units/samples，列为 features 
s <- svd( Y - rowMeans(Y) )
for(i in 1:nrow(Y) ){
 plot(pc$x[,i], s$d[i]*s$v[,i])
}

总结

x <- t(e)
pc <- prcomp(x) # sdev  rotation  center  scale  x
plot(pc$x[, 1], pc$x[, 2], col = group, main = "PCA", xlab = "PC1", ylab = "PC2")

此PCA等效于对centered数据进行SVD，其中居中发生在列（此处为基因）上。

cx <- sweep(x, 2, colMeans(x), "-")
sv <- svd(cx)
plot(sv$u[, 1], sv$u[, 2], col = group, main = "SVD", xlab = "U1", ylab = "U2")

因此，来自SVD的U列对应于PCA中的主要成分x。只是大小不一样。

sv$v[1:5, 1:5]
##            [,1]      [,2]      [,3]      [,4]       [,5]
## [1,]  0.0046966 -0.013275  0.002087  0.017093  0.0006956
## [2,] -0.0021623 -0.002212  0.001543 -0.003346 -0.0034159
## [3,] -0.0030945  0.005870  0.001686  0.003890  0.0019032
## [4,] -0.0007355 -0.002002 -0.002753  0.001776  0.0192205
## [5,]  0.0010133  0.001215 -0.001561  0.003349 -0.0012380
pc$rotation[1:5, 1:5]
##                  PC1       PC2       PC3       PC4        PC5
## 1007_s_at  0.0046966 -0.013275  0.002087  0.017093  0.0006956
## 1053_at   -0.0021623 -0.002212  0.001543 -0.003346 -0.0034159
## 117_at    -0.0030945  0.005870  0.001686  0.003890  0.0019032
## 121_at    -0.0007355 -0.002002 -0.002753  0.001776  0.0192205
## 1255_g_at  0.0010133  0.001215 -0.001561  0.003349 -0.0012380

prcomp返回的standard deviations是sample standard deviations (也就是sample variance的unbiased estimates, 进行了 n/(n−1) correction).
D是主成分的平方和，但不除以sample size.

head(sv$d^2)
## [1] 673418 285393 182527 127667 108576  81999
head(pc$sdev^2)
## [1] 3582.0 1518.0  970.9  679.1  577.5  436.2
head(sv$d^2/(ncol(e) - 1))
## [1] 3582.0 1518.0  970.9  679.1  577.5  436.2

如果没有进行居中，则在以样本为列的数据矩阵上的SVD（与替代变量分析SVA中所使用的）等效于在以样本为行的数据矩阵上的SVD。

sv2 <- svd(t(e))
plot(sv2$u[, 1], sv2$u[, 2], col = group, main = "samples vs genes (typical PCA)", xlab = "U1", ylab = "U2")
sv1 <- svd(e)
plot(sv1$v[, 1], sv1$v[, 2], col = group, main = "genes vs samples (SVA)", xlab = "V1", ylab = "V2")

如果为了比较样本距离，行是样本，基因居中。 为了找到有助于批处理的基因，如在SVA模型中，行是基因，样本进行中心化。