PCA主成分分析实战和可视化 | 附R代码和测试数据

一文看懂PCA主成分分析中介绍了PCA分析的原理和分析的意义(基本简介如下，更多见博客)，今天就用数据来实际操练一下。（注意：用了这么多年的PCA可视化竟然是错的！！！）

在公众号后台回复**“P****CA实战”**，获取测试数据。

一、PCA应用

# 加载需要用到的R包
library(psych)
library(reshape2)
library(ggplot2)
library(factoextra)

1. 数据初始化

# 基因表达数据
exprData <- "ehbio_salmon.DESeq2.normalized.symbol.txt"
# 非必须
sampleFile <- "sampleFile"

2. 数据读入

# 为了保证文章的使用，文末附有数据的新下载链接，以防原链接失效
data <- read.table(exprData, header=T, row.names=NULL,sep="\t")
# 处理重复名字，谨慎处理，先找到名字重复的原因再决定是否需要按一下方式都保留
rownames_data <- make.names(data[,1],unique=T)
data <- data[,-1,drop=F]
rownames(data) <- rownames_data
data <- data[rowSums(data)>0,]
# 去掉方差为0 的行，这些本身没有意义，也妨碍后续运算
data <- data[apply(data, 1, var)!=0,]

3. 数据预处理（可选）

# 计算中值绝对偏差 (MAD, median absolute deviation)度量基因表达变化幅度
# 在基因表达中，尽管某些基因很小的变化会导致重要的生物学意义，
# 但是很小的观察值会引入很大的背景噪音，因此也意义不大。
# 可以选择根据mad值过滤，取top 50， 500等做分析
mads <- apply(data, 1, mad)
data <- data[rev(order(mads)),]
# 此处未选
# data <- data[1:500,]
dim(data)

4. PCA分析

# Pay attention to the format of PCA input
# Rows are samples and columns are variables
data_t <- t(data)
variableL <- ncol(data_t)if(sampleFile != "") {sample <- read.table(sampleFile,header = T, row.names=1,sep="\t")data_t_m <- merge(data_t, sample, by=0)rownames(data_t_m) <- data_t_m$Row.namesdata_t <- data_t_m[,-1]
}
# By default, prcomp will centralized the data using mean.
# Normalize data for PCA by dividing each data by column standard deviation.
# Often, we would normalize data.
# Only when we care about the real number changes other than the trends,
# `scale` can be set to TRUE.
# We will show the differences of scaling and un-scaling effects.
pca <- prcomp(data_t[,1:variableL], scale=T)# sdev: standard deviation of the principle components.
# Square to get variance
# percentVar <- pca$sdev^2 / sum( pca$sdev^2)# To check what's in pca
print(str(pca))

5. PCA结果展示

# PCA结果提取和可视化神器
# http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/
library(factoextra)

1. 碎石图展示每个主成分的贡献

# 如果需要保持，去掉下面第1,3行的注释
#pdf("1.pdf")
fviz_eig(pca, addlabels = TRUE)
#dev.off()

2. PCA样品聚类信息展示

# repel=T，自动调整文本位置
fviz_pca_ind(pca, repel=T)

3. 根据样品分组上色

# 根据分组上色并绘制
fviz_pca_ind(pca, col.ind=data_t$conditions, mean.point=F, addEllipses = T, legend.title="Groups")

4. 增加不同属性的椭圆

“convex”: plot convex hull of a set o points.
“confidence”: plot confidence ellipses around group mean points as the function coord.ellipse() [in FactoMineR].
“t”: assumes a multivariate t-distribution.
“norm”: assumes a multivariate normal distribution.
“euclid”: draws a circle with the radius equal to level, representing the euclidean distance from the center. This ellipse probably won’t appear circular unless coord_fixed() is applied.

# 根据分组上色并绘制95%置信区间
fviz_pca_ind(pca, col.ind=data_t$conditions, mean.point=F, addEllipses = T, legend.title="Groups", ellipse.type="confidence", ellipse.level=0.95)

5. 展示贡献最大的变量 (基因)

1）展示与主坐标轴的相关性大于0.99的变量 (具体数字自己调整)

# Visualize variable with cos2 >= 0.99
fviz_pca_var(pca, select.var = list(cos2 = 0.99), repel=T, col.var = "cos2", geom.var = c("arrow", "text") )

2）展示与主坐标轴最相关的10个变量

# Top 10 active variables with the highest cos2
fviz_pca_var(pca, select.var= list(cos2 = 10), repel=T, col.var = "contrib")

3）展示自己关心的变量（基因）与主坐标轴的相关性分布

# Select by names
# 这里选择的是MAD值最大的几个基因name <- list(name = c("FN1", "DCN", "CEMIP","CCDC80","IGFBP5","COL1A1","GREM1"))
fviz_pca_var(pca, select.var = name)

6. biPLot同时展示样本分组和关键基因

# top 5 contributing individuals and variablefviz_pca_biplot(pca,fill.ind=data_t$conditions,palette="joo",addEllipses = T,ellipse.type="confidence",ellipse.level=0.95,mean.point=F,col.var="contrib",gradient.cols = "RdYlBu",select.var = list(contrib = 10),ggtheme = theme_minimal())

二、PCA分析注意事项

一般说来，在PCA之前原始数据需要中心化（centering，数值减去平均值）。中心化的方法很多，除了平均值中心化（mean-centering）外，还包括其它更稳健的方法，比如中位数中心化等。
除了中心化以外，定标 (Scale, 数值除以标准差) 也是数据前处理中需要考虑的一点。如果数据没有定标，则原始数据中方差大的变量对主成分的贡献会很大。数据的方差与其量级成指数关系，比如一组数据(1,2,3,4)的方差是1.67，而(10,20,30,40)的方差就是167,数据变大10倍，方差放大了100倍。
但是定标(scale)可能会有一些负面效果，因为定标后变量之间的权重就是变得相同。如果我们的变量中有噪音的话，我们就在无形中把噪音和信息的权重变得相同，但PCA本身无法区分信号和噪音。在这样的情形下，我们就不必做定标。
一般而言，对于度量单位不同的指标或是取值范围彼此差异非常大的指标不直接由其协方差矩阵出发进行主成分分析，而应该考虑对数据的标准化。比如度量单位不同，有万人、万吨、万元、亿元，而数据间的差异性也非常大，小则几十大则几万，因此在用协方差矩阵求解主成分时存在协方差矩阵中数据的差异性很大。在后面提取主成分时发现，只提取了一个主成分，而此时并不能将所有的变量都解释到，这就没有真正起到降维的作用。此时就需要对数据进行定标(scale)，这样提取的主成分可以覆盖更多的变量，这就实现主成分分析的最终目的。但是对原始数据进行标准化后更倾向于使得各个指标的作用在主成分分析构成中相等。对于数据取值范围不大或是度量单位相同的指标进行标准化处理后，其主成分分析的结果与仍由协方差矩阵出发求得的结果有较大区别。这是因为对数据标准化的过程实际上就是抹杀原有变量离散程度差异的过程，标准化后方差均为1，而实际上方差是对数据信息的重要概括形式，也就是说，对原始数据进行标准化后抹杀了一部分重要信息，因此才使得标准化后各变量在主成分构成中的作用趋于相等。因此，对同度量或是取值范围在同量级的数据还是直接使用非定标数据求解主成分为宜。
中心化和定标都会受数据中离群值（outliers）或者数据不均匀（比如数据被分为若干个小组）的影响，应该用更稳健的中心化和定标方法。
PCA也存在一些限制，例如它可以很好的解除线性相关，但是对于高阶相关性就没有办法了，对于存在高阶相关性的数据，可以考虑Kernel PCA，通过Kernel函数将非线性相关转为线性相关，关于这点就不展开讨论了。另外，PCA假设数据各主特征是分布在正交方向上，如果在非正交方向上存在几个方差较大的方向，PCA的效果就大打折扣了。

参考资料

https://www.zhihu.com/question/20998460
http://blog.csdn.net/zhongkelee/article/details/44064401
http://www.xiaolingzi.com/?p=963
http://wenku.baidu.com/link?url=hsnzR5gUvsPBwkwwcWU4T3aTSC_fsZDxAmaGGBPfumsIW_I0TJAdJEWhFyiQgw7uA58DKukR-9g5x0DyzE97kHddMaXOxk_iZBjoIdbdB6e
http://stackoverflow.com/questions/22092220/plot-only-y-axis-but-nothing-else
http://www.sthda.com/english/wiki/scatterplot3d-3d-graphics-r-software-and-data-visualization
http://dong.farbox.com/32
http://www.nlpca.org/pca_principal_component_analysis.html
http://gastonsanchez.com/how-to/2014/01/15/Center-data-in-R/
http://www.statpower.net/Content/310/R Stuff/SampleMarkdown.html
http://www2.edu-edu.com.cn/lesson_crs78/self/02198/resource/contents/ch_05/ch_05.html
http://www.sthda.com/english/wiki/principal-component-analysis-in-r-prcomp-vs-princomp-r-software-and-data-mining
http://stackoverflow.com/questions/1249548/side-by-side-plots-with-ggplot2
https://satijalab.org/seurat/pbmc3k_tutorial.html
http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/

PCA主成分分析实战和可视化 | 附R代码和测试数据相关推荐

基于PCA主成分分析的BP神经网络回归预测MATLAB代码
基于PCA主成分分析的BP神经网络回归预测MATLAB代码代码注释清楚. 先对数据集进行主成分分析,自主根据贡献率选择主成分:同时计算KMO验证值:用PCA以后数据进行BP神经网络回归预测. 可以读 ...
独家 | 一文读懂最大似然估计(附R代码)
作者:阿尼·辛格翻译: 陈之炎校对:丁楠雅本文约4200字,建议阅读10+分钟. 本文将研究MLE是如何工作的,以及它如何用于确定具有任何分布的模型的系数. 简介解释模型如何工作是数据科学中最 ...
主成分分析法（PCA）的理解（附python代码案例）
目录一.PCA简介二.举个例子三.计算过程(公式) 3.0 题干假设 3.1 标准化 3.2 计算协方差矩阵 3.3 计算特征值和特征值向量 3.3 多重共线性检验(可跳过) 3.4 适合性检验 ...
python 数据分析可视化实战超全附完整代码数据
代码+数据:https://download.csdn.net/download/qq_38735017/87379914 1.1 数据预处理 1.1.1 异常值检测 ①将支付时间转为标准时间的过程中 ...
python大数据实战项目_商业数据分析比赛实战，内附项目代码
如果你对商业数据分析感兴趣.想要积累更多项目经验,那么就来看看下面这项目吧. 数据竞赛平台和鲸社区最近正在举办一场数据分析大赛,不仅带来了22w奖金和30w创业基金支持,更是提供了统一的在线比赛环境, ...
PCA主成分分析实战案例
遇到的问题: X = df.loc[:,0:4].values#提取第0-3列 y = df.loc[:,4].values #提取第4列报错: TypeError: cannot do slice ...
PCA主成分分析的原理解释及python代码实现
主成分分析(PCA)原理详解 Matlab-PCA Sort eigenvalues and associated eigenvectors after using numpy.linalg.ei ...
180123 PCA主成分分析的原理解释及python代码实现
主成分分析(PCA)原理详解 Matlab-PCA Sort eigenvalues and associated eigenvectors after using numpy.linalg.eig ...
origin三元相图_三元相图怎么看怎么画（附R代码示例）
什么是三元相图三元图是重心图的一种,它有三个变量,但需要三者总和为恒定值.在一个等边三角形坐标系中,图中某一点的位置代表三个变量间的比例关系.在群体遗传学中,它被称做Finetti图:在博弈论中,常 ...