主成分分析（Principal component analysis，PCA）是将多个存在相关关系的变量，降维为几个互不相关的指标，也就是主成分。
我的数据相对来说，有点复杂，188个品种，11个表型性状，性状间的相关性也很强，直接去做后续分析的话，数据比较庞大，也不太好解释，所以在这里引入主成分分析，来简化数据。
读取数据
第一列为品种，第2到12列为11个性状。
> dat <- read.table("blue.txt", header = T, check.names = F, sep = "\t")
> head(dat)
Cul UFW NR AL UDW TL PA SA AD V
1 1 0.3243333 22.03333 3.310000 0.03300000 46.85667 3.214000 10.098000 0.6920667 0.1807000
2 2 0.3177726 15.55898 2.538101 0.03366908 32.98066 2.778663 8.727610 0.8714627 0.1892379
3 3 0.3335468 18.45504 2.400459 0.03783424 37.59174 2.986453 9.310346 0.8395869 0.1894390
4 4 0.2413333 20.23333 1.955000 0.02800000 24.11700 2.124000 6.674333 0.8038667 0.1542667
5 5 0.3090000 19.60000 1.950000 0.03633333 32.47000 2.474667 7.767000 0.8270000 0.1557000
6 6 0.3336667 23.00000 1.550333 0.03766667 27.44200 2.604333 8.175000 0.9568333 0.1995667
NT MTL
1 54.00000 65.02167
2 30.76526 44.28370
3 40.44400 47.85501
4 35.00000 35.87500
5 36.40000 47.51000
6 31.66667 38.41167
1. FactoMineR包PCA（）函数
1.1 主成分分析
PCA分析的模型如下：
PCA(data, scale.unit = TRUE, ncp = 5 , graph = TRUE)
# scale.unit：布尔值（即只有TRUE和FALSE两种），默认为TRUE，会在进行分析之前对数据进行标准化
# ncp：数字，最终结果中保留的维数，默认为5。这个ncp我没有太懂，如果理解为只保留5个维度，后面的PCA还是会保留和变量数目一样的维数。
# graph：布尔值，如果为TURE，则会在分析完成之后直接输出PCA的图，但由于我们会在图上进行修改和美化，因此本例中使用FALSE
提取dat中第2到12列表型测量值进行PCA，并输出PCA结果，其中包括了以下组件，用来描述PCA的各项指标。其中，var是variables的缩写，ind是individuals的缩写，后面的内容是对前面指标的解释。
> library('FactoMineR')
> library("ggplot2")
> library("factoextra")
> pca <- PCA(dat[,2:12],graph = FALSE)
> print(pca)
**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 188 individuals, described by 11 variables
*The results are available in the following objects:
name description
1 "$eig" "eigenvalues"
2 "$var" "results for the variables"
3 "$var$coord" "coord. for the variables"
4 "$var$cor" "correlations variables - dimensions"
5 "$var$cos2" "cos2 for the variables"
6 "$var$contrib" "contributions of the variables"
7 "$ind" "results for the individuals"
8 "$ind$coord" "coord. for the individuals"
9 "$ind$cos2" "cos2 for the individuals"
10 "$ind$contrib" "contributions of the individuals"
11 "$call" "summary statistics"
12 "$call$centre" "mean of the variables"
13 "$call$ecart.type" "standard error of the variables"
14 "$call$row.w" "weights for the individuals"
15 "$call$col.w" "weights for the variables"
1.2 提取主成分特征值和累计贡献率
> eig.val <- get_eigenvalue(pca)
> eig.val
eigenvalue variance.percent cumulative.variance.percent
Dim.1 7.121720585 64.74291441 64.74291
Dim.2 1.843421579 16.75837799 81.50129
Dim.3 1.075634876 9.77849887 91.27979
Dim.4 0.428591056 3.89628233 95.17607
Dim.5 0.166737398 1.51579452 96.69187
Dim.6 0.126993851 1.15448955 97.84636
Dim.7 0.108295415 0.98450377 98.83086
Dim.8 0.081041291 0.73673901 99.56760
Dim.9 0.035208509 0.32007735 99.88768
Dim.10 0.010527211 0.09570192 99.98338
Dim.11 0.001828229 0.01662026 100.00000
eigenvalue：第i个主成分的特征值
variance.percent：第i个主成分的贡献率
cumulative.variance.percent：第1到第i个主成分的累积贡献率
1.4 确定主成分个数
指标1：通常选择特征值四舍五入，大于1的主成分，本数据中，选取前三个。
指标2：累积贡献率，通常累积贡献率达到80%以上即可，在本数据中，前两个主成分的累积贡献率就已经达到了81.50%。
指标3：碎石图。根据每个主成分的贡献率（variance.percent）绘制柱形图，当贡献率降低的趋势基本稳定，降低到一定程度时，选择拐点之前的主成分个数。如下图所示，本数据中应该选取三个主成分。
综合以上三个指标，我最终确定选择三个主成分。
> fviz_eig(pca, addlabels = TRUE,ylim = c(0, 65))
1.4 可视化主成分分析结果
#获得主成分分析中各变量的分析结果
> var <- get_pca_var(pca)
> var
Principal Component Analysis Results for variables
===================================================
Name Description
1 "$coord" "Coordinates for the variables" #变量做散点图的坐标
2 "$cor" "Correlations between variables and dimensions" #变量和主成分之间的相关性
3 "$cos2" "Cos2 for the variables" #变量在因子图中代表的变异的质量
4 "$contrib" "contributions of the variables" #变量对主成分的贡献率
主成分1中，各变量的贡献率：
从下图中可以明显看出，主成分1中，比较重要的变量分别为SA、PA、TL、MTL、NT、UFW和V。
> fviz_contrib(pca, choice ="var", axes = 1)
同理，主成分2中，做出主要贡献的变量是AD和AL。
> fviz_contrib(pca, choice ="var", axes = 2)
主成分3的主要变量是AL和NR。
> fviz_contrib(pca, choice ="var", axes = 3)
各变量在主成分1、2、3中总体的贡献量，下图中，红线为各变量贡献量的均值，超过均值则被认为在三个主成分中具有重要贡献。
> fviz_contrib(pca, choice ="var", axes = 1:2:3)
1.5 变量坐标(coord)与相关性(cor)可视化
根据分组信息绘制：
> set.seed(123) #设定三个主成分
> res.km <- kmeans(var$coord, centers = 3, nstart = 25)
> grp <- as.factor(res.km$cluster)
> fviz_pca_var(res.pca, col.var = grp,
+ palette = c("#0073C2FF", "#EFC000FF", "#868686FF"),
+ legend.title = "Cluster")
2. psych包principal()函数
FactoMineR的优势在于可视化美观方便，但是好像不能直接输出主成分得分，因此选择psych包中的principal()函数进行主成分得分计算，并和FactoMineR的结果进行比较。
principal(r, nfactors = 1, residuals = FALSE,
rotate="varimax",n.obs=NA, covar=FALSE,
scores=TRUE,missing=FALSE,
impute="median",oblique.scores=TRUE,
method="regression",...)
r：指定输入的数据，如果输入的是原始数据，R将自动计算其相关系数矩阵；
nfactors：指定主成分个数；
residuals：是否显示主成分模型的残差，默认不显示；
rotate：指定模型旋转的方法，默认为最大方差法；
n.obs：如果输入的数据是相关系数矩阵，则必须指定观测样本量
covar：为逻辑参数，如果输入数据为原始数据或方阵（如协方差阵），R将其转为相关系数矩阵；
scores：是否计算主成分得分；
missing：缺失值处理方式，如果scores为TRUE，且missing也为TRUE，缺失值将被中位数或均值替代；
impute：指定缺失值的替代方式，默认为中位数替代；
method：指定主成分得分的计算方法，默认使用回归方法计算。
2.1 确定主成分个数
> library(psych)
> fa.parallel(x = dat, fa="pc")
看图应该选择2个主成分，但是结合实际数据，还有贡献率以及特征值，这里确定为三个主成分。
2.2 提取主成分
> pc <- principal(r = dat[,2:12], nfactors = 3, rotate = 'none')
# rotate = 'none'：不进行主成分旋转
> pc
Principal Components Analysis
Call: principal(r = dat[, 2:12], nfactors = 3, rotate = "none")
Standardized loadings (pattern matrix) based upon correlation matrix
PC1 PC2 PC3 h2 u2 com
UFW 0.87 0.37 0.11 0.89 0.106 1.4
NR 0.77 0.23 -0.54 0.94 0.063 2.0
AL 0.30 -0.62 0.69 0.95 0.047 2.4
UDW 0.69 0.27 0.36 0.67 0.326 1.9
TL 0.95 -0.29 -0.05 0.98 0.019 1.2
PA 0.99 0.05 0.02 0.97 0.026 1.0
SA 0.99 0.06 0.02 0.98 0.023 1.0
AD -0.09 0.89 0.34 0.92 0.078 1.3
V 0.86 0.35 0.11 0.88 0.118 1.4
NT 0.88 -0.28 -0.16 0.89 0.113 1.3
MTL 0.92 -0.33 -0.04 0.96 0.040 1.3
# 各个变量在三个主成分中的载荷量，根据系数，可以写出各主成分得分的计算公式。
# h2为主成分对每个变量的方差解释度
# u2为主成分无法解释的比例
PC1 PC2 PC3
SS loadings 7.12 1.84 1.08 #特征值
Proportion Var 0.65 0.17 0.10
Cumulative Var 0.65 0.82 0.91
Proportion Explained 0.71 0.18 0.11 #贡献率
Cumulative Proportion 0.71 0.89 1.00 #累积贡献率
Mean item complexity = 1.5
Test of the hypothesis that 3 components are sufficient.
The root mean square of the residuals (RMSR) is 0.03
with the empirical chi square 16.45 with prob < 0.9
Fit based upon off diagonal values = 1
可以看到结果和FactoMineR基本一致，其中特征值完全相同，但是贡献率FactoMineR较高。
根据载荷量可以写出公式：
X1=0.87*UFW+0.77*NR+0.30*AL+0.69*UDW+0.95*TL+0.99*PA+0.99*SA-0.09*AD+0.86*V+0.88*NT+0.92*MTL
X2=0.37*UFW+0.23*NR-0.62*AL+0.27*UDW-0.29*TL+0.05*PA+0.06*SA+0.89*AD+0.35*V-0.28*NT-0.33*MTL
X3=0.11*UFW-0.54*NR+0.69*AL+0.36*UDW-0.05*TL+0.02*PA+0.02*SA+0.34*AD+0.11*V-0.16*NT-0.04*MTL
2.3 计算各主成分的得分
> spc <- principal(r = dat[,2:12], nfactors = 3, rotate = 'varimax', scores = TRUE)
> scores <- data.frame(spc$scores)
> head(scores)
RC1 RC2 RC3
1 0.7545669 -0.1792164 0.3835327
2 0.1627193 1.6261451 0.0191181
3 0.4799638 1.3190902 -0.2408160
4 -0.2400794 0.7343836 -1.1764561
5 0.1905985 1.0733563 -0.8608670
6 0.3709758 2.1503158 -1.5103191
> write.table(scores, file = "pca_scores.txt", row.names = F, col.names = T, quote = F, sep = "\t")
2.4 确定各主成分权重
各主成分在最终主成分得分里所占的权重是不同的，通过贡献率来计算：
Wi=variance.percent(i) /cumulative.variance.percent(n)
n：最终主成分确定个数
即第i个主成分的权重=第i个主成分的贡献率/第i到第n个主成分的累积贡献率
参考资料：
PCA - Principal Component Analysis Essentials - Articles - STHDA
https://mp.weixin.qq.com/s/oNKqEjJtiJYsAR6ALrPVrQ
R语言中的PCA分析与可视化_nikang3148的博客-CSDN博客
以上内容为个人粗浅的理解，欢迎讨论，如有错误，也敬请指出。
网友评论