PCA的核心思想：假设我们的数据有D1, D2,···,Dn个维度，PCA就是要构造一个线性变换，PCi = W1D1 + W2D2 +···+ WnDn，Wj就是第j维度【列】在第i个【行】PC中的权重。找PC【主成分】有先后顺序，我们总是先找总方差最大的PC，方差解释度是统计里最重要的一个概念。然后我们再找与前一个PC线性无关的能解释最多方差的下一个PC，以此类推，知道得到所有的n个PC。最终的结果就是原先的n个维度通过线性变换，变成了新的n个线性无关的按方差解释度排序的PC。主成分中的“主principal”针对的就是方差解释度。
PCA的输入很简单，就是一个矩阵，行为样本，列为维度/特征；
数学 语文 英语
小明
小红 Wj
小方
小强
这里我把scale这一步单独列出来了，大家可以看下scale的效果是什么，PCA的原始数据如果没有scale，PCA的结果就是错的，因为PCA非常看重方差占比，我们的数据的单位往往千差万别，必须scale，使得不同维度的方差可比。
PCA的输出，sdev、rotation、x；sdev标准差是每个PCA的方差解释度的平方根，肯定是递减的；rotation旋转，很形象，PCA只是把原维度做了旋转而已，rotation是个矩阵里面就是原维度如何通过线性变换，得到PC的；x就是PCA之后原来数据在新的PC下的坐标。
my_data <- mtcars[,c(1,3)]
head(my_data)
mpg disp
Mazda RX4 21.0 160
Mazda RX4 Wag 21.0 160
Datsun 710 22.8 108
Hornet 4 Drive 21.4 258
Hornet Sportabout 18.7 360
Valiant 18.1 225
plot(my_data, main=“raw data”)
my_data <- scale(my_data, center = T, scale = T)
plot(my_data, main = “scaled data”)
######### mtcars.pca <- prcomp(mtcars[,c(1:7,10,11)], center = TRUE,scale. = TRUE)
mtcars.pca <- prcomp(my_data, center = F,scale. = F)
summary(mtcars.pca)
Importance of components:
PC1 PC2
Standard deviation 1.3592 0.39045
Proportion of Variance 0.9238 0.07622
Cumulative Proportion 0.9238 1.00000
plot(mtcars.pca$x, main = “after PCA”)
sdev
the standard deviations of the principal components (i.e., the square roots of the eigenvalues of the covariance/correlation matrix, though the calculation is actually done with the singular values of the data matrix).
rotation
the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.
x
if retx is true the value of the rotated data (the centred (and scaled if requested) data multiplied by the rotation matrix) is returned. Hence, cov(x) is the diagonal matrix diag(sdev^2).
For the formula method, napredict() is applied to handle the treatment of values omitted by the na.action.
输入输出明确了，接下来我们开始做一些个性化的事情：
- 如果我有新的数据，如何看它在原PCA里的位置？
看下PC是如何得到的？？？？，从rotation里我们可以得到线性变换的方程。
PC1 PC2
mpg -0.7071068 0.7071068
disp 0.7071068 0.7071068
PC1 = -0.7071068 * mpg + 0.7071068 * disp
PC2 = 0.7071068 * mpg + 0.7071068 * disp
把新数据代入到这个方程里就可以求新数据的PC score了，在R里，一般用矩阵乘法实现：Newdata %*% rotation
- 如何找对我们某个PC贡献最大或top10的维度？取绝对权重最大的top维度集合
在对RNA-seq数据进行PCA分析时非常实用，我们发现了感兴趣的PC，往往要找哪些基因贡献最大，或者直接可视化我们某个基因在前两个PC中的位置。
蓝线就代表了我们原坐标下的PC1方向，如何反向找出蓝线的斜率呢？
假设蓝线为y=kx，那么向量[1, k]在rotation的作用下必然变为[c, 0]，c的大小和符号是不确定的，但我们有一个确定的0，有点线性代数基础的就知道如何解k了。
plot(my_data, main = “scaled data”)
abline(v=0)
abline(h=0)
abline(coef = c(0, -0.7071068/0.7071068), col=“blue”)
abline(coef = c(0, 0.7071068/0.7071068), col=“red”)
那么我们就知道了第一个方差最大的方向是y = -0.7071068/0.7071068 * x，第二个方差最大的就是y = 0.7071068/0.7071068.
普通实现
如何用普通方法，不用线性代数和矩阵来求呢？
先明确总方差这个概念：
单一维度的方差好理解，直接求var就行，比如我们scale后，mpg和disp的方差就都为1了。
但是我们的数据PCA之后，PC1的方差就为1.847425，PC2的方差为0.1524486。
Importance of components:
PC1 PC2
Standard deviation 1.3592 0.39045
Proportion of Variance 0.9238 0.07622
Cumulative Proportion 0.9238 1.00000
如何理解总方差呢？其实很简单，就是所有PC的方差之和。PC1的Proportion of Variance是这么算的：1.847425/(1.847425 + 0.1524486)
现在我们手动来做PCA，先求第一个PC，就是要找蓝线y = kx，如何求k？
我们要把所有点投射到蓝线上，然后要让所有的投影点的在蓝线这个维度方差最大化。
这里没办法了，点太多不可能手算，所以只能用最大似然估计了。
最大似然估计用的mle函数，此函数需要一个似然函数，变得就是A，就是Ax+By+C=0，我们就是不停地试A，然后看什么情况下我们所有的点在此线上投影的方差最大。
最后的结果是：1.009254*x + y = 0，可以看到与原来的线几乎重合，我们手动算出来的PC和prcomp算出来的一模一样。
这里的这个似然函数的构造有一点技巧：需要知道怎么算点到线的距离；根据勾股定理算投影长度，需要判断投影的正负。
PC1定了，其实PC2也就定了，因为只有一条线与PC1正交，就是垂直的线。所以也就不重复算了。
PCA的线性代数实现
自此，我们已经了解了PCA的思想内涵，普通实现。
这个时候再来讲PCA的线性代数实现才是合适的，为什么线性代数里的工具能高效的实现PCA？
R的PCA包，prcomp和princomp，The function princomp() uses the spectral decomposition approach. The functions prcomp() and PCA()[FactoMineR] use the singular value decomposition (SVD).
先从直觉理解PCA。先看一个数据实例，明显的两个维度之间有一个相关性，大部分的方差可以被斜对角的维度解释。
我们的任务就是通过线性变化，找出一个维度能最大的解释数据的方差；然后继续找正交的第二个维度，直至最后一个正交维度。
不用矩阵思维我们都可以做，构造一个线性组合，优化方案就是方差最大，然后寻找另一个正交维度。我们的目的是将covariacne matrix对角化，也就是变量内的方差最大，变量间的相关程度协方差为0；特征向量的矩阵表征的就是线性变化，所以原先的数据矩阵在施加x运动后，会投射到新的PC里，这也从直觉上理解了PCA的精髓，线性投射。
从定义出发，Ax=cx：A为矩阵，c为特征值，x为特征向量。矩阵A乘以x表示，对向量x进行一次转换（旋转或拉伸）（是一种线性转换），而该转换的效果为常数c乘以向量x（即只进行拉伸）。我们通常求特征值和特征向量即为求出该矩阵能使哪些向量（当然是特征向量）只发生拉伸，使其发生拉伸的程度如何（特征值大小）。这样做的意义在于，看清一个矩阵在那些方面能产生最大的效果（power），并根据所产生的每个特征向量（一般研究特征值最大的那几个）进行分类讨论与研究。知乎Rex
talia有个PPT讲得非常严谨。
不得不惊叹数学的神奇，一个复杂问题用这么精简巧妙的方式解决了。
PCA，以基因表达数据为例，每个基因在样本间都是在变得，而且有很多是相互相关的，举个例子，如果两个/一群基因是完全相关的，那我们完全只需要一个变量/维度就可以解释样本里的所有变异。PCA就是通过线性变化，把原来的所有变量（n个）转化为新的n个PC，这些PC彼此间线性无关，而且可以根据变异的解释量来排序PC，每个变量对新的PC都有一个简单的线性转换公式。
例子：以小学生基本生理属性为案例，分享下R语言的具体实现，分别选取身高（x1）、体重(x2)、胸围(x3)和坐高(x4)。具体如下：
student<- data.frame(
- x1=c(148,139,160,149,159,142,153,150,151),
- x2=c(41 ,34 , 49 ,36 ,45 ,31 ,43 ,43, 42),
- x3=c(72 ,71 , 77 ,67 ,80 ,66 ,76 ,77,77),
- x4=c(78 ,76 , 86 ,79 ,86 ,76 ,83 ,79 ,80)
- )
student
x1 x2 x3 x4
1 148 41 72 78
2 139 34 71 76
3 160 49 77 86
4 149 36 67 79
5 159 45 80 86
6 142 31 66 76
7 153 43 76 83
8 150 43 77 79
9 151 42 77 80
student.pr <- princomp(student,cor=TRUE)
student.pr
Call:
princomp(x = student, cor = TRUE)
Standard deviations:
Comp.1 Comp.2 Comp.3 Comp.4
1.8846028 0.5738007 0.3094410 0.1525488
4 variables and 9 observations.
summary(student.pr,loadings=TRUE)
Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4
Standard deviation 1.884603 0.57380073 0.30944099 0.152548760
Proportion of Variance 0.887932 0.08231182 0.02393843 0.005817781
Cumulative Proportion 0.887932 0.97024379 0.99418222 1.000000000
Loadings:
Comp.1 Comp.2 Comp.3 Comp.4
x1 0.510 0.436 0.139 0.728
x2 0.513 -0.172 0.741 -0.398
x3 0.473 -0.761 -0.396 0.201
x4 0.504 0.448 -0.524 -0.520
screeplot(student.pr,type=“lines”)
standard deviation：标准差
Proportion of Variance：方差的占比
Cumulative Proportion：累计贡献率
由上图可见四项指标做分析后，给出了4个成分，他们的重要性【方差的占比】分别为：0.887932、0.08231182、0.02393843、0.005817781，累积贡献【累计贡献率】为：0.887932、0.97024379、0.99418222 1.000000000各个成分的碎石图也如上，可见成份1和成份2的累积贡献已经达到95%，因此采用这两个成份便可充分解释学生的基本信息。
可以算出 Z1 和 Z2的公式
temp<-predict(student.pr)
plot(temp[,1:2])
主成分分析（Principal Component Analysis，PCA）的目标是用一组较少的不相关的变量代替大量相关变量，同时尽可能保留原始变量的信息，推导所得的变量就成为主成分，是原始变量的线性组合。也就是将N个变量（N维），通过线性组合，降维到K个综合变量（K维，K <N）来归纳性解释某一个现象。
先举个简单例子帮助理解吧：
某篮球俱乐部有40名男同学，同学之间各个指标存在或大或小的差异，包括身高、体重、视力、百米速度、肺活量、每天练球时长、睡眠时间等指标。在季度选拔中，40名同学进球得分数（成绩）存在差异，那么这些指标是否均与成绩相关？或者相关性有多大？
此时可能就要用到PCA，分析方法简述如下：
1选择初始变量
比如以上7个指标作为变量（a1 - a7），40名同学作为样本。
2对原始数据矩阵进行标准化，做相关系数矩阵
（1）原始数据矩阵：每行为40名男同学的各项指标值，每列为各项指标在40名同学中的体现；
（2）因为各个指标度量单位不同，取值范围不同，不宜直接由协方差矩阵出发，因此选择相关系数矩阵。
3计算特征值及相应特征向量
4判断主成分的个数
最常见的方式是根据特征值判断，一般选择特征值大于1的变量数作为PC个数，假如，此项分析中特征值大于1的有两个，则最终可以有2个主成分（具体主成分个数可以根据实际研究调整）。
5得到主成分表达式
通过分析得到第一主成分（PC1）和第二主成分（PC2），假设表达式如：（a1*表示a1标准化后的数值）
PC1=0.92a1+0.89a2+…
pc2=…
6结合数据的实际意义开展分析
PC1比PC2更能解释样本间差异的原因（如下图中横纵轴的百分比）。PC1的线性组合中a1、a2、a3、a6贡献度较大（前面的系数较大），PC2的线性组合中a5贡献度较大。以PC1为横轴，各个样本根据成绩大小有明显区分，说明以上7个指标中，身高、体重、视力、每天练球时长这4个指标与同学的成绩相关性更强。为了迎合线性组合的概念，应该找个更合适的词语来综合描述身高、体重、视力、每天练球时间这4个指标，以涵盖这4个指标的意义。理解了以上的概念，再将主成分分析用于RNA-seq也很容易理解。
扫一扫
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
