R包做主成分分析（含实例）

注：内容非原创，来自对网站http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/的翻译。

1 什么是主成分分析？

Principal component analysis is used to extract the important information from a multivariate data table and to express this information as a set of few new variables called principal components. 

主成分分析（Principal component analysis，PCA）属于多元统计分析中的一种，适用于将具有一定相关性的多个变量化为少数几个相互无关的综合指标（主成分），以方便可视化。

比如我这里有两个地区采集到的多个A虫样本（individuals），测量了31个形态指标(即变量variables)，我想知道两个地区的A虫是否存在形态上的差异，就用我的数据做一下主成分分析。

真正理解主成分分析原理需要线性代数的知识，这里只作简单解释：

如下图A，数据本来在xy坐标系中显示出来，通过找到数据展现出最大变异的PC1，再把数据投影到PC1上，我们实现了降维，把在二维坐标系中的数据降到了一维。PC2轴是第二重要的方向，它和PC1正交，也就是说二者之间不相关。

The amount of variance retained by each principal component is measured by the so-called eigenvalue.

每个主成分保留的变异又通过特征值来表示。

在数据集中的变量彼此高度相关时，主成分分析尤其有用。变量相关就说明数据存在冗余。

上图来感受下什么是高度冗余（两个变量高度相关）：

什么是较低的冗余（两个变量相关不明显）：

2 R中如何实现主成分分析

我们通过R包"FactoMineR"和"factoextra"来实现主成分分析和可视化。

下面我会通过自己的一个实例来展示。

#安装和载入所需要的R包

>install.packages(c("FactoMineR","factoextra","ggplot2"))

>library("FactoMineR")

>library("ggplot2")

>library("factoextra")

“FactoMineR”是用来实现数据分析的，“factoextra”是用来实现基于“ggplot2”的可视化的

首先展示一下我的数据：

前两列是样本编号（样本，individuals）和每个样本的采样地点（sites），不参与数据分析。后面31列对应的是31个变量（观测值或者变量，variables）

#从剪切板读取数据，前提是在excel里面已经复制数据集

>collina=read.table("clipboard",header=TRUE)

#显示数据框的前几行，确保正确读取数据

>head(collina)

#进行主成成分分析，第1、2列不参与数据分析（第一列是巢编号，第二列是采样地点）

>collina.pca <- PCA(collina[,-c(1,2)],graph = FALSE)

一般在数据分析之前要进行标准化处理（standardization），尤其是当变量的测量单位不同的时候（如千米，厘米等）。标准化的目的是为了让变量可比，总的来说就是标准差为1，平均值为0。

进行标准化时，数据可进行如下转化：

用R自带的scale()可以进行标准化，但FactoMineR里的PCA()已经包括了标准化的步骤，因此不需要额外的标准化处理

3 对主成分分析的解读和可视化

3.1 主成分的特征值（eigenvalues）和方差（variance）

前面说过，特征值衡量的是主成分保留的变异量，第一主成分保留的变异量最大，依次递减。第一主成分也与数据集方差最大的方向相同。通过特征值，我们可以决定主成分的个数。

#提取主成分（PCs）的特征值和保留方差的百分比（信息），输入：

>eig.val <- get_eigenvalue(collina.pca)

>eig.val

出来的结果如下图：

第一列是每个主成分的特征值。

第二列是每个特征值解释的变异的百分比，计算方法是某主成分特征值/特征值的加和。

第三列是累计的变异的百分比

在本例中，第一、二主成分能解释48.9%的变异。

我们在前面说过特征值可以用来决定保留的主成分的个数，那么怎么来确定呢？

方法1 选择特征值大于1 的主成分保留

在数据已经进行标准化处理的前提下，主成分特征值大于1，意味着主成分比标准化后的原始数据集的任一变量解释的方差要多。根据这个标准，本例中可以保留前7个主成分。

方法二 根据累计方差百分比

一般认为主成分的累计方差百分比达到80%以上，就可以了。根据这个标准，本例中可以保留前7个主成分。

当然80%这个值是可变的，如果你认为70%就够了也是可以的。

方法3 做陡坡图来确定

陡坡图就是将主成分按照第二列其，特征值解释的变异的百分比的大小顺序来排列的柱状图。当趋势降低到已经变得差不多小，没有太大减小的趋势时，认为前面的主成分就是可以保留的。

#做陡坡图，确定主成分的个数

>fviz_eig(collina.pca, addlabels = TRUE,ylim = c(0, 32))

结果如下图：

根据陡坡图，本例中，可以保留前4个主成分。

3.2获得和可视化变量结果

3.2.1 获得主成分分析中变量的结果

#获得主成分分析中变量的结果

>var <- get_pca_var(collina.pca)

>var

这个命令可以帮助我们得到关于变量的主成分分析的结果，帮助我们进行变量的作图，具体帮助如下：

 var$coord:coordinates of variables to create a scatter plot

var$coord提供用变量做散点图的坐标

var$cos2: represents the quality ofrepresentation for variables on the factor map. It’s calculated as the squaredcoordinates: var.cos2 = var.coord * var.coord.

var$cos2代表的是变量在因子图中代表的变异的质量

var$contrib: contains the contributions(in percentage) of the variables to the principal components. The contributionof a variable (var) to a given principal component is (in percentage) :(var.cos2 * 100) / (total cos2 of the component).

var$contrib代表的是变量对于主成分的贡献（百分比表示）。

#查看前4个变量对主成分的贡献

>head(var$contrib, 4)

# 变量对主成分1（PC1）的贡献作图，显示贡献前20的变量

fviz_contrib(collina.pca, choice ="var", axes = 1, top = 20)

# 变量对主成分1（PC1）的贡献作图，显示贡献前20的变量

fviz_contrib(collina.pca, choice ="var", axes = 2, top = 20)

# 变量对主成分1、2（PC1和PC2）的贡献作图，显示贡献前26的变量

fviz_contrib(res.pca, choice ="var", axes = 1:2, top = 26)

其中红色的虚线代表的是期望的平均贡献，如果变量的超过了这条线，认为是对主成分有重要贡献的。

3.2.2 获得主成分分析中样本（individuals）的结果

>ind <- get_pca_ind(res.pca)

>ind

具体含义和变量的结果类似。

要查看各项结果输入：

# Coordinates of individuals

head(ind$coord)

# Quality of individuals

head(ind$cos2)

# Contributions of individuals

head(ind$contrib)

3.2.3 biplot作图

The representation of variables differs from the plot of the observations: The observations are represented by their projections, but the variables are represented by their correlations.

主成分分析的结果主要通过biplot图来呈现，横坐标是主成分1，纵坐标是主成分2，变量通过与主成分的相关关系表现在图中，样本点通过该样本多个变量的观测值在主成分上的投影表现在图中。输入：

>fviz_pca_biplot(collina.pca,

                #individuals

                geom.ind = "point",#只显示样本点，不包括文本信息，即样本编号

                fill.ind = collina$sites,#按sites进行分组

                pointshape = 21, pointsize = "contrib",#点的尺寸大小表示样本的贡献

                palette = "jco",#配色选择“jco”杂志的调色板

                addEllipses = TRUE,#在相同组的样本点为画出密度椭圆

                #variables

                col.var = "cos2",#变量的颜色表示其在因子图中代表的变异的质量，即cos2

                gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),repel=TRUE,

                legend.title = list(fill = "sites"))

出来的PCA分析的biplot图如下：

我们可以看出来，变量的cos2值的大小通过不同的颜色梯度表示出来，每个样本点对主成分的贡献通过点的大小表现出来。两个采样地的样本点也进行了分组。

但是，由于我自身数据的原因，两个采样点的样本分不太开。