ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
122 stars 30 forks source link

PCA-Statistics is the new sexy!!! #576

Closed ixxmu closed 3 years ago

ixxmu commented 3 years ago

https://mp.weixin.qq.com/s/lSoEBXYMCuxp93Vi7DEvBw

github-actions[bot] commented 3 years ago

PCA-Statistics is the new sexy!!! by 生信技能树

标题看着熟悉的,都是看过《神探夏洛克》的,真真儿是爱这张马脸;
  • 一切随缘吧,感觉这篇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:231: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

免费的课程和学习资源 :