|
最近我们分享了不少降维的分析，例如PCoA，NMDS，CCA等。前几天有同学在公众号留言，想学习一下主成分分析(PCA)。这一期我与大家分享下PCA分析的基本知识，以及如何利用R语言计算PCA的相关参数，最后我把数据处理及其可视化的整个过程打包成函数，便于新手直接套用，同时也希望能给已经入门的同学在学习自定义函数的过程提供一定的帮助。 咱们常提到的PCA分析全称是Principal Components Analysis，即主成分分析，这是降维中最常见的一种方法。我相信，数据维度大家应该都比较清楚，这里就不再一一赘述。首先，我们要清楚明白什么是降维？直白的说就是把数据的维度降下来，用相对低维的向量来表征高维度数据的特征。其次，为什么要做降维呢，降维有什么必要性以及好处使得我们需要适用数据更低的维度表示？大背景是因为大数据时代，多维的数据越来越普遍，比如组学和微生物OTU的数据，数据维度都成千上万，处理起来太复杂。通过降维的方法找出最具代表性的主成分，便于后续对样本进行区分和数据处理。 数据降维的必要性：
1）降低高维数据的处理难度，便于后续计算和分析；（也就是说，100个维度通过降维形成几个主要维度，并用其来表征100个维度做后续的分析） 2）去除噪音和冗余数据，减少主要信息的损失；（维度降下来意味着只保留最主要的规律和信息，那些细小的相关性是噪音的影响） 3）保证计算的准确性和效率；（维度越高相对来说的确精度越低，一是数据度量本身的不准确性增加，二是计算时比如浮点数或者舍入等情况越多的发生） 另外，PCA是一种无监督算法，意味着我们不需要标签就能对数据进行降维；降维后，由于失去了标签，我们可能无法理解每个维度的含义；但至少减少了数据维度，使得计算机能更好的识别和计算。也就是说，PCA把原来高维的数据（多个特征）用低维（少量特征值）来代替，新的维度是原先高维度的线性组合，这些组合变量（方差最大）尽可能代表原来的变量，而且彼此之间互不相关，因此对于一些冗余的数据有很好的表现。（如果大家想了解降维的一些原理，推荐大家阅读https://zhuanlan.zhihu.com/p/74501834） 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分析才是咱们今天的重头戏~ #逐步计算主要参数
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函数的用法简单，与常规的求取特征值和特征向量不同的是，prcomp函数是对变量矩阵(相关矩阵)采用SVD方法计算其奇异值（原理上是特征值的平方根）。
例子如下： #相关矩阵分解
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 #所有样本每个轴的得分
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分析，我利用基础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))结果如下图： 
|
