PCA-Statistics is the new sexy!!!

标题看着熟悉的，都是看过《神探夏洛克》的，真真儿是爱这张马脸；

  • 一切随缘吧，感觉这篇PCA，大家不一定看得进去；前天的脚本命名里说到佛度有缘人，那就看看有没有同样感兴趣的有缘人吧。

  • 把基于基础函数的PCA放在文章里了，能够加深理解，希望你喜欢。

了解PCA

PCA的重要性，我昨天已经介绍了：PCA-弱水三千，取哪一瓢饮？

PCA是为了更好地展示多维数据，通过线性转化，展示保留最多信息的主成分；将样本尽可能地分散地展示在坐标轴中达到可视化的目的；

  • PCA的理论假设是：方差越大，信息量越大；

  • 拿生信数据来说，大概率上，我们是要看数据的分组情况，所以要在坐标轴上表示的point是样本，舍掉的是基因层面的component，即n个基因获得n个component，依据方差最大化，取前k(0<k<n)个component;

  • 本质上计算出n个特征向量，给予矩阵n个移动方向，最后保存了k个移动后的结果；

PCA步骤：

1）数据为m行n列的原始矩阵（sample为行，gene为列）
2）矩阵X每一个元素减去该列的均值（中心化）
目的是使所有维度的偏移都是以0为基础的（我们必须对数据中individual（sample）和observations（gene）有区分和了解）
3）求出协方差矩阵
4）目的是协方差矩阵中除对角线外的元素为0，即实现协方差矩阵对角化；
5）将P按特征值进行排序，因为Y=PX，所以，中心化后的矩阵（转置）与特征向量矩阵（转置）乘积后得到新的样本矩阵，取前两行即PC1和PC2；

这里把PCA的过程用我理解的基础函数，做了包装，大家试着理解一下吧：

rm(list=ls())
######数据集可用于测试PCA
library("FactoMineR")
library("factoextra")
data("decathlon2")
decathlon2.active <- decathlon2[1:23, 1:10]
decathlon2.active

pca_base<-function(data,x=1,y=2){
  center_d<-scale(data,center=TRUE,scale=FALSE)
  cov_deca<-cov(center_d)
  deca_rotation<-eigen(cov_deca)
  PC<- (t(deca_rotation$vectors)%*%t(center_d))[x:y,]
  return(PC)
}
pca_base(data = decathlon2.active)

我们汲汲以求的PCA其实早有对统计学烂熟于心的人做了R包，不得不说，数学才是王道啊！！！

对比下在R的现成的PCA功能的结果

  • FactoMineR和factoextra配合做PCA和可视化(下图中图片名为PCA)；

  • prcomp（stats base级别）和autoplot配合做PCA和可视化(下图中图片名为prcomp)；

######以下是FactoMineR和factoextra的工作：
res<-PCA(X = decathlon2.active, scale.unit = FALSE, ncp = 5, graph = FALSE)
res$ind$coord###PCA图中采用的坐标
######fviz_pca_ind对individual绘图，fviz_pca_var对variable绘图
fviz_pca_ind(res)

######prcomp为stats自带的PCA方法
deca_re<-prcomp(decathlon2.active)
#####rotation-包含特征向量的矩阵
deca_re$rotation[, 1]
#####x-如果retx参数设为TRUE，则返回rotated矩阵（中心化（scaled，如果有设定）矩阵乘以rotation矩阵）的结果；cov（x）为协方差矩阵；
deca_re$x
###生成名为Rplot_deca$x的图，方便与ggfortify的作图对比
plot(deca_re$x)
#####ggfortify使ggplot2功能更加丰富，使autoplot能够处理prcomp的结果
library(ggfortify)
autoplot(prcomp(decathlon2.active),label=TRUE,loadings=TRUE)

########两个PCA方法对比
#####对coord处理后获得特征向量，与prcomp中的rotation一致
loadings<-sweep(res$var$coord,2,sqrt(res$eig[1:5,1]),FUN="/")
loadings

PCA
看着跟前面的图坐标位置哪儿哪儿不一样，后面再用$x画图后，看到是坐标比例的差异，再对比发现，与上图是某种镜像的关系，相对位置其实是一样的：
prcomp
Rplot_deca$x

参考内容：
1.https://www.zhihu.com/question/36481348
2.http://www.sthda.com/english/wiki/print.php?id=204
3.https://wenku.baidu.com/view/ce7ee04bcc175527072208ca.html
4.https://zhuanlan.zhihu.com/p/21580949
5.https://groups.google.com/forum/#!topic/factominer-users/BRN8jRm-_EM
6.http://blog.sciencenet.cn/home.php?mod=space&uid=443073&do=blog&id=321347

赞 (0)

相关推荐