主成分分析笔记

主成分分析简介及主要思想

主成分分析（principal components analysis，PCA），是由1901年被Pearson首先引入的，1933年由Hotelling作了进一步的发展，主成分分析是从多个数值变量（指标）之间的相互关系入手，利用降维的思想，将多个变量（指标）化为少数几个互不相关的综合变量（指标）的统计方法。

在医学研究中，为了客观、全面地分析问题，常要记录多个观察指标并考虑众多的影响因素，这样的数据虽然可以提供丰富的信息，但同时也使用数据的分析工作更趋复杂化。例如在儿童生长发育的评价中，收集到的数据包括每一个儿童的身高、体重、胸围、坐高、肺活量等十多个指标。怎样利用这类多指标的数据对每一个儿童的生长发育水平作出正常的评价？如果仅用其中任一指标来作评价，其结论显然是片面的，而且不能充分利用已有的数据信息。如果分别利用每一指标进行评价，然后再综合各指标评价的结论，这样做一是可能会出现各指标评价的结论不一致，甚至相互冲突，从而给最后的综合评价带来困难；二是工作量明显增大，不利于进一步的统计分析。事实上，在实际工作中，所涉及的众多指标之间经常是有相互联系和影响的，从这一点出发，希望通过对原始相互关系的研究，找出少数几个综合指标，这些综合指标是原始指标的纯属绵且，它既保留了原始指标的主要信息，且又互不相关。这样一种从众多原始指标之间相互关系入手，寻找少数综合指标以根据原始指标信息的多元统计方法称为主成分分析。

主成分分析的数学模型及几何意义

（一）主成分分析的数学模型

设有m个指标，即$X_{1},X_{2},\dots,X_{m}$,欲寻找可以概括这m个指标主要信息的综合指标$Z_{1},Z+{2},\dots,Z_{m}$，从数学上讲，就是寻找一组常数，即$a_{i1},a_{i2},\dots,a_{im}(i=1,2,\dots,m）$，使这m个指标的线性组合为（以下为公式1-1——：

$$\left\{\begin{array}{c}{Z_{1}=a_{11} X_{1}+a_{12} X_{2}+\cdots+a_{1 m} X_{m}} \\ {Z_{2}=a_{21} X_{1}+a_{2} X_{2}+\cdots+a_{2 m} X_{m}} \\ {\vdots} \\ {Z_{m}=a_{m 1} X_{1}+a_{m 2} X_{2}+\cdots+a_{m m} X_{m}}\end{array}\right.$$

能够根据m个原始指标$X_{1},X_{2},\dots,X_{m}$的主要信息，且各$Z_{i}(i=1,2,\dots,m)$互不相关，为叙述方便，可以引入以下的矩阵形式，令

则公式1-1可以表示为：

$$\bf {Z=AX} (1-2)$$

或者是

如果$Z_{1}=\bf a’_{1}X$满足$\bf a’_{1}a_{1}=1$，且
$Var(Z_{1})=\mathop{Max}\limits_{\bf a’a=1}^{}{Var(\bf a’X)}$，刚称$Z_{1}$是原始指标$X_{1},X_{2},\dots,X_{m}$的第一主成分。

通常情况下，如果$Z_{i}=\bf a’X$满足：
(1)$\bf a’_{i}a_{i}=1$，当i>1时，$\bf a’_{i}a_{j}=0(i=1,2,\dots,i-1)$
(2)$Var(Z_{i})=\mathop{Max{Var(a’\bf X)}}\limits_{\bf a’a=1,\bf a’a_{j}=0(j=1,2,\dots,i-1)}^{}$
刚称$Z_{i}$是原始指标的第$i$主成分（$i=2,\dots,m$）。

由上述定义可知，当$i\neq j$时，主成分$Z_{i}$与$Z_{j}$是互不相关的，并且$Z_{1}$是原始指标$X_{1},X_{2},\dots,X_{m}$的一切线性组合中方差最大者，$Z_{2}$是与$Z_{1}$不相关的、除$Z_{1}$以外的$X_{1},X_{2},\dots,X_{m}$一切线性组合中方差最大者，$Z_{m}$是与$Z_{1},Z_{2},\dots,Z_{m}-1$都不相关的，除$Z_{1},Z_{2},\dots,Z_{m}-1$以外的$X_{1},X_{2},\dots,X_{m}$一切线性组合中方差最大者。从理论上讲，求得的主成分个数最多可能有m个，这时，m个主成分就反映了全部原始指标所提供的信息，鉴于主成分分析的目的主要是用较少的综合指标来反映全部原始指标中的主要信息，因此在实际工作中，所确定的主成分个数总是小于原始指标的个数。

（二）主成分的几何意义

为方便讨论，以m=2为例说明主成分分析的几何意义，设个体具有两个观测指标X1和X2，它们之间具有较强的相关性，测量n例这样的个体的值，将所得的n对数据在以X1为横轴，X2为纵轴的二维坐标平面中的苫，得到如下的散点图：

由上图可以看出，由于$X_{1}$和$X_{2}$具有较强的相关性，这n个点的分页呈现出直接化的趋势；同时它们沿$X_{1}$轴方向和$X_{2}$轴方向都具有较大的变异变。个体在某个方向上的变异度可以用该方向上相应观测变量的方差来定量地表示。显然，如果只考虑$X_{1}$和$X_{2}$中任何一个方向上的方差，就将损失原始观测数据中很大一部分信息。如果我们将坐标轴$X_{1}$和$X_{2}$同时按逆时针方向作一个放置，得到新的坐标轴$Z_{1}$和$Z_{2}$，使得在亲折坐标平面上，这n个点的分布基本上不再具有有相关性，且它们的变异主要集中在Z1方向上，而在Z2方向上的变异较小，此时若取$Z_{1}$作为第一主成分，则$Z_{1}$就反应了原始指标$X_{1}$和$X_{2}$所包含的主要信息。

主成分的求法及性质

（一）主成分的求法

由主成分的定义可知，各主成分互不相关，即任意两个主成分$Z_{i},Z_{j}$的协方差为0，即

$$Cov(Z_{i},Z_{j})=0,i\neq j$$

且各主成分的方差满足：

$$Var(Z_{1}) \geq Var(Z_{2}) \dots Var(Z_{m})$$

于是由公式(1-2)定义的随机向量$Z$的协方差矩阵为：

由主成分定义中的条件（1）可知，这里的方阵$A$是正交阵，即$A’A=I$（I为单位矩阵），由此可得：

由上述公式可知，求原始指标$X_{1},X_{2},\dots,X_{m}$的主成分问题，实际上就是要求满足上述条件的正交阵$A$，即随机微量$X=(X_{1},X_{2},\dots,X_{m}$的协方差矩阵$Cov(X)$的特征值（eigenvalue）与特征向量(eigenvector)。

主成分的计算过程

下面讨论怎么由一组$X_{1},X_{2},\dots,X_{m}$的样本观测值求出主成分，假设收集到的原始数据共有n例，每例测得m个指标的数值，记录如下所示：

1. 对各原始指标数据进行标准化，先按下式进行：

$$X'_{ij}=\frac{X_{ij}-\bar{X_{j}}}{S_{J}}(j=1,2,3,\dots,m)$$

将原始指标标准化（标准化通俗讲就是每一列的每个数字每去这一列的均值，然后除以标准差），然后用标准化的数据$X’_{ij}$来计算主成分，为了方便计算，下面的公式中仍用$X_{ij}$表示标准化后的指标数据，$\bf X$为标准化后的数据矩阵，则：

2. 求出$X$的协方差矩阵$R$

标准化后，$X$的相关矩阵即为协方差矩阵$Cov({\bf X})$

对角线上分别是$r_{11}$，$r_{22}$和$r_{mm}$的方差，非对角线上是协方差。协方差是衡量两个变量同时变化的变化程度。以两个变量x、y来举例，例如协方差大于0，表示x和y若一个增，另一个也增；小于0表示一个增，一个减。如果ｘ和ｙ是统计独立的，那么二者之间的协方差就是０；但是协方差是０，并不能说明ｘ和ｙ是独立的。协方差绝对值越大，两者对彼此的影响越大，反之越小。协方差是没有单位的量，因此，如果同样的两个变量所采用的量纲发生变化，它们的协方差也会产生数字上的变化。

3. 求出相关矩阵的特征值和特征值对应的特征向量

由公式1-3得知，求主成分的问题，实际上是求出$X$的协方差矩阵$Cov(X)$（这里即为$X$的相关矩阵$R$）的特征值和特征向量，由于$R$为半正定矩阵，故可由R的特征方程

$$|{\bf R} -\lambda{\bf I|} = 0$$

解得每一特征值$\lambda_{i}$对应的单位特征向量$a_{i}=(a_{i1} a_{i2} \dots a_{im}’$，从而求得各主成分，即

$$Z_{i}={\bf a'_{i}X}={a_{i1}X_{1}+a_{i2}X_{2}+ \dots +a_{im}X_{m}} i=1,2,\dots,m$$

有关特征值和特征向量的理解可以参考后文引用的文章。这里简单提一下，如果把矩阵看作是运动，对于运动而言，最重要的就是运动的速度和方向，那么特征值就是运动的速度，特征向量就是运动的方向。

（二）主成分的性质

1. 各主成分互不相关

即$Z_{i}$与$Z_{j}$的相关系数为0，即：

$$r_{z_{i},z_{j}}=\frac{Cov(Z_{i},Z_{j})}{\sqrt{Cov(Z_{i},Z_{j})Cov(Z_{j},Z{j})}}=0(i \neq j)$$

因此各主成分间的相关系统矩阵为单位矩阵。

2. 主成分的贡献率和累积贡献率

可以证明，各原始指标$X_{1},X_{2},\dots,X_{m}$的方差和与各主成分$Z_{1},Z_{2},\dots,Z_{m}$的方差和相等，将数据标准化后，原始指标的方差和为$\sum\limits_{i=1}^m \lambda_{i}$，即有$m=\sum\limits_{i=1}^m \lambda_{i}$

各指标所提供的信息量是用其方差来衡量的。由此可知，主成分分析是把m个原始指标$X_{1},X_{2},\dots,X_{m}$的总方差分解为m个互不相关的综合指标$Z_{1},Z_{2},\dots,Z_{m}$的方差之和，使第一主成分的方差达到最大，即变化最大的方向微量所相应的线性函数，最大方差为$\lambda$。其中${\lambda}/{\sum\limits_{i=1}^m \lambda_{i}}$表明了第一主成分$Z_{1}$的方差在全部方差中所占的比值，称为第一主成分的贡献率，这个值越大，表明$Z_{1}$这个指标综合原指标$X_{1},X_{2},\dots,X_{m}$的能力越强，也可以说，由$Z_{1}$的差异来解释$X_{1},X_{2},\dots,X_{m}$的差异的能力越强，正是因为这一点，才把$Z_{1}$称为$X_{1},X_{2},\dots,X_{m}$的第一主成分，也就是$X_{1},X_{2},\dots,X_{m}$的主要部分，了解到这一点，就可以知道为什么主成分是按特征值$\lambda_{1},\lambda_{2},\dots,\lambda_{m}$进行排序的。

通常情况下

$$\frac{\lambda_{i}}{\sum\limits_{i=1}^m}=\frac{\lambda_{i}}{m}(k=1,2,\dots,m)$$

为第i主成分的贡献率；而称

$$\sum\limits_{i=1}^k \frac{\lambda_{i}}{m}(k \leq m)$$

为前k个主成分的累积贡献率。

3. 主成分个数的选取

通常并不需要全部的主成分，只用其中的前几个，一般来说，主成分的保留个数按以下的原则来进行：

(1)以累积贡献率来确定：

当前k个主成分的累积贡献率达到某一特定值时（一般以大于70%为宜，有的时候会要求大于80%），则保留前k个主成分。

(2)以特征值的大小来确定：

即若主成分$Z_{i}$的特征值$\lambda_{i} \geq 1$，则保留$Z_{i}$，否则就去掉该主成分。这个与碎石图类似。

4. 因子载荷

为了了解各主成分与各原始指标之间的关系，在主成的表达式中，第$i$主成分$Z_{i}$的特征值的平方根$\sqrt{\lambda_{i}}$与第$j$原始指标$X_{j}$的系数$a_{ij}$的乘积，即

$$q_{ij}=\sqrt{\lambda_{i}a_{ij}}$$

为因子载荷（factor loading），由因子载荷构成的矩阵称为因子载荷阵，事实上因子载荷$q_{ij}$就是第$i$主成分$Z_{i}$与第$j$原始指标$X_{j}$之间的相关系数，它反映了主成分$Z_{i}$与原始指标X_{i}$之间联系的密切程度与作用的方向。

5. 样品的主成得分

对于具有原始指标测定值$(X_{i1},X_{i2},\dots,X_{im})$的任一样品，可先用标准化变换式将原始数据标化，即：

$$X'_{ij}=\frac{X_{ij}-\bar{X_{j}}}{S_{J}}(j=1,2,3,\dots,m)$$

然后代入各主成分的表达式，即

$$Z_{i}={\bf a'_{i}X}={a_{i1}X_{1}+a_{i2}X_{2}+ \dots +a_{im}X_{m}} i=1,2,\dots,m$$

求出该样本各主成分值，这样求得的主成分值称为该样本的主成分得分，利用样品的主成分得分，可以对样品的特征进行推断和评价。

案例分析

案例一

例22-1 某研究者测得84名10岁男孩的身高、坐高、体重、胸围、肩宽、肺活量等6项生长发育指标，数据见表22-2。试作主成分分析。

第1步：导入数据

raw_pca <- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_szq_2201.csv")
# 这是原始数据，放在了github上
head(raw_pca)
str(raw_pca)
raw_pca <- raw_pca[-1]
#第1列是序号，没用，剔除掉

剔除掉第1列后的数据：

第2步： 数据的简单统计量

apply(raw_pca,2,sd) # 求各列的标准差
apply(raw_pca,2,mean) # 求各列的均数

第3步：相关矩阵的特征值

raw.pr<- prcomp(raw_pca, scale = TRUE) # 标准误、方差贡献率和累积贡献率
summary(raw.pr,loadings=TRUE) # 标准误、方差贡献率和累积贡献率

在上图中，其中第二行，即Proportion of Variance表示的是贡献率；而第3行，即Cumulative Proportion表示的是累积贡献率；由这个表可以看出，主成分取3个比较合适，此时的累积贡献率为88.92%，接近于90.00%。

第4步：相关矩阵的特征向量

prcomp(raw_pca, scale = TRUE)

这张表的结果可以得出各个主成分的与相应变量的系数，如果取3个主成分，则可以得出前三个主成分为：

第5步：因子载荷矩阵

前面的第1到第5步用的是R中的base的函数，没有找到计算因载荷的矩阵方法，现在用psych包来进行计算。

碎石图

# 利用paych包进行分析
library(psych)
fa.parallel(raw_pca,fa='pc',n.iter=100,show.legend = FALSE,main='PCA analysis')

psych包中的fa.parallel用来提取主成分，绘制的图形如下所示：

这种图叫碎石图，把对应各个主成分的特征值按从大到小的顺序在图上绘出，选取主成分个数至碎石图发生斜率明显变化为止。其中，蓝线是基于观测特征值的碎石检验，根据100个随机数据矩阵推导出来的特征值均值是红线，从图中可以看出：碎石检验图形最大变化处上面只有一个成分；特征值大于随机模拟数据的也只有一个主成分；特征值大于1的也只有一个主成分。所以这一组数据使用一个主成分即可保留数据集的大部分信息（这一点与书上的不符，但主成分的挑选我觉得参入了一定的主观成分在里面，选取几个主成分，要综合几个因素来看，书上选取的是3个主成分，为了保持一致，我们后面的分析也选3个主成分）。

提取主成分

提取3个主成分，如下所示：

raw2.pca <- principal(r = raw_pca, nfactors = 3, rotate = 'none')

# nfactors=3表示，提取3个主成分
# rotate指定旋转的方法，默认是最大方差旋转；scores设定是否需要计算主成分得分，默认不计算。
raw2.pca

结果如下所示：

红色框中就是因子载荷矩阵。

由因子载荷矩阵可知，第一主成分$Z_{1}$在各原始指标上的因子载荷较为均匀，故可认为该主成分反映的是各原始指标的综合信息；第二主成分$Z_{2}$在$X_{1}$（身高）、$X_{3}$（坐高）及$X_{4}$（胸围）上的因子载荷较大，故可认为该主成分反映的是体型方面的信息；而第三主成分$Z_{3}$则主要反映了来自原始指标$X_{6}$（肺活量）的信息。此外，还可知道，第一主成分$Z_{1}$与各原始变量之间的关系较为密切，第二主成分$Z_{2}$与原始变量$X_{1}$、$X_{2}$及$X_{4}$之间的关系较为密切，而第三主成分$Z_{3}$与原始变量$X_{6}$之间的关系较为密切。

此外，结果中还有h2和u2以及com，它们的解释如下：
h2栏指成分公因子方差，主成分对每个变量的方差解释度；
u2栏指成分唯一性，就是方差不能被主成分解释的比例。显然u2=1-h2；

主成分旋转

当我们提取的主成分不止一个时，使用主成分旋转会使各主成分所代表的实际意义更容易被解释。主要的旋转方法包括以下两种：

1. 正交旋转：选择的成分保持不相关
2. 斜交旋转：选择的成分保持相关
   在本案例中，我们选择使用正交旋转中的方差极大旋转，它试图对载荷阵的列进行去噪，使得每个成分只由一组有限的变量来解释（即载荷阵每列只有少数几个很大的载荷，其他都是很小的载荷）。

# varimax:方差极大旋转
raw_pca_xz <- principal(raw_pca,nfactors = 3,rotate = 'varimax')
raw_pca_xz

经过主成分旋转的三个主成分仍然不相关，对变量的解释性不变，这是因为变量的群组没有发 生变化。另外，三个主成分旋转后的累积方差解释性没有变化（89%），变的只是各个主成分对方差的解释度（第一主成分从71%变为39%，第二主成分从10%变为30%，第三主成分从8%变为20%）。

获取主成分与原始变量的线性关系（保存在模型的weights部分）：

head(raw_pca_xz$weights)

这样，我们就可以得到主成分与原始变量之间的相关关系：

RC1 = -0.18X1 - 0.273X2 + 0.47X3 + …
RC2 =0.535X1 + 0859X2 - 0.067X3 + …
RC3 = 0.049X1 - 0.267X2 - 0.173X3 + …

代码总结与美化

上面的代码只是为了说明PCA的原理，在实际绘图中可以使用其他的R包来进行PCA图的生成，下面是一段生成PCA的代码：

## README
# 0. Set working directory 
rm(list=ls())
temp_path <- file.path("D:","FIGURE","PCA")
dir.create(temp_path,recursive = TRUE)
# dir.create() function need to add recursive = TRUE when create multi-level directory

setwd(temp_path)
getwd()


# 1. Generation test dataset
test_dataset <- matrix(nrow=100, ncol=10)

data.matrix <- test_dataset
colnames(data.matrix) <- c(
  paste("wt",1:5,sep=""),
  paste("ko",1:5,sep=""))
rownames(data.matrix) <- paste("gene", 1:100, sep="")
for (i in 1:100){
  wt.values <- rpois(5, lambda = sample(x=10:1000, size=1))
  ko.values <- rpois(5, lambda = sample(x=10:1000, size=1))
  
  data.matrix[i,] <- c(wt.values, ko.values)
}

head(data.matrix)

# If you have data for PCA analysis 
# Plaease Transform your data into type same to data.matrix


# 2. Calculate PCA
pca <- prcomp(t(data.matrix),scale=TRUE)
plot(pca)
pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100,1)

# 3. Calcualte loading scores

loading_score <- pca$rotation[,1]
# prcomp() calcuate results inclding loading score, namely, rotation term

gene_scores <- abs(loading_score)
# Some numbers in loading _score is negative
# Use abs() function transform them into positive

gene_score_ranked <- sort(gene_scores, decreasing = TRUE)
# Sort all loading score after abs()


top_10_genes <- names(gene_score_ranked[1:10])
top_10_genes
# Obtain top10 genes in loading score

pca$rotation[top_10_genes,1]
# Examine top10 gene of loading score

# 3. Plot each component  for Scree plot

barplot(pca.var.per, 
        main="Scree Plot", 
        xlab="Principal Component", 
        ylab="Percent Variation")

# 4. ggplot2 foir PCA image
pca.data <- data.frame(Sample=rownames(pca$x),
                       X=pca$x[,1],
                       Y=pca$x[,2])
# one column with the sample ids
# Tow columns for the X and Y coordinates for each sample
pca.data

library(ggplot2)
ggplot(data=pca.data, aes(x=X, y=Y, label=Sample))+
  geom_text() +
  xlab(paste("PC1 - ", pca.var.per[1], "%", sep=""))+
  ylab(paste("PC2 - ", pca.var.per[2], "%", sep=""))+
  theme_bw()+
  ggtitle("My PCA Graph")

# 5. Save reulst
ggsave("PCA.pdf",
       device = "pdf",
       dpi = 300,
       limitsize = FALSE
       )

参考资料

1. 医学统计学.孙振球.第4版
2. 如何通俗易懂地解释「协方差」与「相关系数」的概念？
3. 浅谈协方差矩阵
4. 如何理解矩阵特征值和特征向量？
5. 数据挖掘系列篇(15):R语言VS主成分分析的案例
6. Learn R | 数据降维之主成分分析（下）