【R分享|实战】PCA分析与可视化
主成分分析(PCA)的基本介绍
咱们常提到的PCA分析全称是Principal Components Analysis，即主成分分析，这是降维中最常见的一种方法。我相信，数据维度大家应该都比较清楚，这里就不再一一赘述。首先，我们要清楚明白什么是降维？直白的说就是把数据的维度降下来，用相对低维的向量来表征高维度数据的特征。其次，为什么要做降维呢，降维有什么必要性以及好处使得我们需要适用数据更低的维度表示？大背景是因为大数据时代，多维的数据越来越普遍，比如组学和微生物OTU的数据，数据维度都成千上万，处理起来太复杂。通过降维的方法找出最具代表性的主成分，便于后续对样本进行区分和数据处理。
数据降维的必要性：
1）降低高维数据的处理难度，便于后续计算和分析；（也就是说，100个维度通过降维形成几个主要维度，并用其来表征100个维度做后续的分析）
2）去除噪音和冗余数据，减少主要信息的损失；（维度降下来意味着只保留最主要的规律和信息，那些细小的相关性是噪音的影响）
3）保证计算的准确性和效率；（维度越高相对来说的确精度越低，一是数据度量本身的不准确性增加，二是计算时比如浮点数或者舍入等情况越多的发生）
另外，PCA是一种无监督算法，意味着我们不需要标签就能对数据进行降维；降维后，由于失去了标签，我们可能无法理解每个维度的含义；但至少减少了数据维度，使得计算机能更好的识别和计算。也就是说，PCA把原来高维的数据（多个特征）用低维（少量特征值）来代替，新的维度是原先高维度的线性组合，这些组合变量（方差最大）尽可能代表原来的变量，而且彼此之间互不相关，因此对于一些冗余的数据有很好的表现。（如果大家想了解降维的一些原理，推荐大家阅读https://zhuanlan.zhihu.com/p/74501834）
PCA常用的参数
PCA分析中需要用到几个常用参数，标准化(scale)、特征值(eigen value)、特征向量(eigen vector)、载荷(loading)、得分(score)：
1）标准化
如果是针对环境因子，各变量之间存在不同量纲，标准化可以较好地解决这个问题；其次，有些数据的数值比较大，标准化后可以较好地避免较大的数值对主成分的贡献过大。尽管如此，标准化也不是十全十美的。标准化会导致各变量之间的权重相等，这可能会产生一些负面结果。特别是，如果变量中有噪音的话，标准化无形中把噪音和信息的权重变得相同，但PCA是无法区分噪音和信号。因此，在这种情形下，我们就不需要标准化了。
2）特征值和特征向量
特征值表示标量部分，一般为某个主成分的方差，其相对比例可理解为方差解释度或贡献度 ；特征值从第一主成分逐渐减小。特征向量为对应主成分的线性转换向量（线性回归系数），特征向量与原始矩阵的矩阵积为主成分的得分。特征向量是单位向量，其平方和为1。
由于比较抽象，借助该文理解，https://zhuanlan.zhihu.com/p/165382601。“从线性代数的角度出发，如果把矩阵看作n维空间下的一个线性变换，这个变换有很多的变换方向，我们通过特征值分解得到的前N个特征向量，就对应了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向，就可以近似这个矩阵（变换）。其中的N个变化方向，就是这个矩阵最重要的“特征”。
清楚特征的概念后，进一步如何理解特征值与特征向量？如果把矩阵看作是位移，那么特征值 = 位移的速度，特征向量 = 位移的方向。
特征向量在一个矩阵的作用下作伸缩运动，伸缩的幅度由特征值确定。特征值大于1时，所有属于此特征值的特征向量变长；特征值属于(0, 1)，特征向量缩短；特征值小于0，特征向量则反向延长。”
3）载荷
因子载荷矩阵并不是主成分的特征向量，即不是主成分的系数。主成分系数的求法：各自因子载荷向量除以各自因子特征值的算数平方根。
4）得分
是指主成分的得分，其求法：矩阵与特征向量的积。
由于本人非数学和统计学专业，关于PCA分析中一些专业词汇的解释没办法做到位。尽我所能来给大家简化和解读了，为了让大家有一个更准确和更深刻的理解，附上一些比较专业的解读希望能够给大家提供一定的帮助。
当然，最主要的还是讲R语言的部分。如何利用R语言完成PCA分析才是咱们今天的重头戏~
R语言如何PCA分析
#逐步计算主要参数
library(tidyverse)
#使用自带数据集"iris"
#特征分解
eigen <- scale(iris[,-5])%>% # scale()对数据进行标准化
cor()%>% # cor()对矩阵进行相关性分析
eigen() # eign()对数据矩阵进行光谱分解，得到特征值和特征向量
eigen
#提取特征值
eigen$values
#特征向量提取
eigen$vectors
#计算标准化的主成分得分
scale(iris[,-5])%*%eigen$vectors%>%
head()
结果如下：
接下来，我们介绍两个被广泛用于PCA分析的函数，分别为prcomp函数和princomp函数~
?prcomp例子如下：
#相关矩阵分解
iris.pca <- prcomp(iris[,-5],
scale. = T, #scale.=T表示标准化
rank. = 4, #rank.指定最大秩的数字
retx=T) #retx一个逻辑值，指示返回已旋转的变量
iris.pca
#查看结果
summary(iris.pca) #方差解释度
iris.pca$sdev #主成分的标准偏差
iris.pca$rotation #可变载荷矩阵
iris.pca$x #所有样本每个轴的得分
3）princomp函数
princomp以计算相关矩阵或者协方差矩阵的特征值为主。
?princomp# prcomp()的用法
iris.pc <- princomp(iris[,-5],cor=T,scores = T)
iris.pc
summary(iris.pc) #各主成份的解释量
iris.pc$loading #可变载荷矩阵
head(iris.pc$score) #所有样本各轴的得分
screeplot(iris.pc,type="lines") #方差分布图
head(predict(iris.pc)) #预测
biplot(iris.pc,scale=F) #直接把x与rotation绘图，而非标准化
由碎石图可以看出，第三个主成分之后，图线变化趋于平稳，因此可以选择前三个主成分做分析（由于这里只有4列变量，所以效果并不明显）。另外，根据上面主成分解释量的结果来看，其实选前两轴数据进行分析即可（因为轴一解释了73%，轴2解释了22.9%）。
根据轴1和轴2的数据绘制了PCA的降维分析图。
自定义函数用于PCA分析及其可视化
为了更便捷且更完美地执行PCA分析，我利用基础prcomp函数以及R包ggbiplot中的ggbiplot函数构建了一个新函数 (定义名为：simple_PCA)。其他内容就不再赘述，具体看代码吧~
# 自定义函数
# 设置了两个参数，对象数据及主成分数量
simple_PCA <- function(data, components = c(1:2)) { #函数主体
library(ggbiplot) #加载ggbiplot包使用ggbiplot函数
# 重新排序分组信息
data[, 1] <- factor(data[, 1], levels = unique(data[, 1]))
# 除去物种信息名，只保留数值型数据进一步做PCA分析
data_PCA <- data[c(2:ncol(data))]
# 利用prcomp函数进行PCA分析
iris.pca <- prcomp(data_PCA,
center = TRUE, # 中心化和标准化
scale. = TRUE)
# 这里大家可以help一下ggbiplot函数查阅一下对应参数用法
# 具体参数过多，这里就不一一细讲，有问题的可以在推文留言
ggbiplot(iris.pca, choices = components,
obs.scale = 1, var.scale = 1,
ellipse = F, groups = data[, 1], # 注意这里即可，刚刚是定义了
#数据集第一列的分组信息，这里保持一致即可，其他参数可以直接使用
ellipse.prob = 0.95, circle = F,
varname.size = 5, varname.adjust = 1.5,
var.axes = T) +
geom_point(aes(shape = groups, colour = groups), size = 3) +
geom_hline(yintercept = 0, colour = "black", linetype = "longdash", size = 1) +
geom_vline(xintercept = 0, colour = "black", linetype = "longdash", size = 1) +
scale_color_discrete(name = '') +
theme(legend.direction = "horizontal", legend.position = "top") +
stat_ellipse(type = "t", size = 1, geom = "polygon", alpha = 0.2, aes(fill = groups), level = 0.85) +
stat_ellipse(type = "t", size = 1, aes(colour = groups), level = 0.85) +
guides(groups = FALSE) +
theme(text = element_text(size = 15)) +
ggtitle("PCA") +
theme(plot.title = element_text(hjust = 0.5))
}
# 由于上述定义的函数中对应数据集的第一列为分组信息
# 这里构建一个以第一列为分组信息的新数据集
data <- data.frame(groups=iris$Species, iris[1:4])
# 查看数据结构
str(data)
# 以轴1和轴2绘图
simple_PCA(data,components = c(1:2))
结果如下图：