PCA原理及其R实现

点击下方图片查看HappyChart专业绘图软件

主成分分析法

主成分分析也称主分量分析，是揭示大样本、多变量数据或样本之间内在关系的一种方法，旨在利用降维的思想，把多指标转化为少数几个综合指标，降低观测空间的维数，以获取最主要的信息。
在统计学中，主成分分析（principal components analysis, PCA）是一种简化数据集的技术。它是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上，依次类推。主成分分析经常用减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。但是，这也不是一定的，要视具体应用而定。

主成分分析法的基本原理

主成分分析法是一种降维的统计方法，它借助于一个正交变换，将其分量相关的原随机向量转化成其分量不相关的新随机向量，这在代数上表现为将原随机向量的协方差阵变换成对角形阵，在几何上表现为将原坐标系变换成新的正交坐标系，使之指向样本点散布最开的p 个正交方向，然后对多维变量系统进行降维处理，使之能以一个较高的精度转换成低维变量系统，再通过构造适当的价值函数，进一步把低维系统转化成一维系统。
#主成分的一般定义
设有随机变量$X_1，X_2，\dots ，X_p$，样本标准差记为$S_1，S_2，\dots ，S_p$。首先作标准化变换：
$$C_j=a_{j1}{x}_1+a_{j2}{x}_2+\dots +a_{jp}{x}_p,j=1,2,\dots ,p$$
我们有如下的定义：
(1) 若$C_1=a_{11}{x}_1+a_{12}{x}_2+\dots +a_{1p}{x}_p$， ，且使 $Var(C_1)$最大，则称$C_1$为第一主成分；

(2) 若$C_2=a_{21}{x}_1+a_{22}{x}_2+\dots +a_{2p}{x}_p，(a_{21}，a_{22}，\dots ，a_{2p})$垂直于$(a_{11}，a_{12}，\dots ，a_{1p})$，且使$Var(C_2)$最大，则称$C_2$为第二主成分；

(3) 类似地，可有第三、四、五…主成分，至多有p个。
#主成分的性质
主成分$C_1，C_2，\dots ，C_p$具有如下几个性质：
(1) 主成分间互不相关，即对任意i和j，$C_i$ 和$C_j$的相关系数
$$Corr(C_i,C_j)=0$$

(2) 组合系数$(a_{i1}，a_{i2}，\dots ，a_{ip})$构成的向量为单位向量，

(3) 各主成分的方差是依次递减的， 即 $$Var(C_1)\ge Var(C_2)\ge \dots \ge Var(C_p)$$
(4) 总方差不增不减， 即
$$Var(C_1)+Var(C_2)+\dots +Var(C_p)=Var(x_1)+Var(x_2)+\dots +Var(x_p)=p$$
这一性质说明，主成分是原变量的线性组合，是对原变量信息的一种改组，主成分不增加总信息量，也不减少总信息量。

(5) 主成分和原变量的相关系数 $Corr(C_i，x_j)=a_{ij} $

(6) 令$X_1，X_2，\dots ，X_p$的相关矩阵为R, $(a_{i1}，a_{i2}，\dots ，a_{ip})$则是相关矩阵R的第$i$个特征向量(eigenvector)。而且，特征值$l_i$就是第i主成分的方差， 即 $Var(C_i)=l_i$
其中$l_i$为相关矩阵R的第i个特征值(eigenvalue)
$$l_1\ge l_2\ge \dots \ge l_p\ge 0$$

主成分数目的选取

前已指出，设有p个随机变量，便有p个主成分。由于总方差不增不减，$C_1，C_2$等前几个综合变量的方差较大,而$C_p，C_{p-1}$等后几个综合变量的方差较小, 严格说来，只有前几个综合变量才称得上主(要)成份，后几个综合变量实为“次 ”(要)成份。实践中总是保留前几个，忽略后几个。
保留多少个主成分取决于保留部分的累积方差在方差总和中所占百分比(即累计贡献率)，它标志着前几个主成分概括信息之多寡。实践中，粗略规定一个百分比便可决定保留几个主成分；如果多留一个主成分，累积方差增加无几，便不再多留。
#主成分分析的主要作用
1．主成分分析能降低所研究的数据空间的维数。即用研究m维的Y空间代替p维的X空间(m＜p)，而低维的Y空间代替 高维的x空间所损失的信息很少。即：使只有一个主成分$Y_l$(即 m＝1)时，这个$Y_l$仍是使用全部X变量(p个)得到的。例如要计算$Y_l$的均值也得使用全部x的均值。在所选的前m个主成分中，如果某个$X_i$的系数全部近似于零的话，就可以把这个$X_i$删除，这也是一种删除多余变量的方法。
2．有时可通过因子负荷aij的结论，弄清X变量间的某些关系。
3．多维数据的一种图形表示方法。我们知道当维数大于3时便不能画出几何图形，多元统计研究的问题大都多于3个变量。要把研究的问题用图形表示出来是不可能的。然而，经过主成分分析后，我们可以选取前两个主成分或其中某两个主成分，根据主成分的得分，画出n个样品在二维平面上的分布况，由图形可直观地看出各样品在主分量中的地位，进而还可以对样本进行分类处理，可以由图形发现远离大多数样本点的离群点。
4．由主成分分析法构造回归模型。即把各主成分作为新自变量代替原来自变量x做回归分析。
5．用主成分分析筛选回归变量。回归变量的选择有着重的实际意义，为了使模型本身易于做结构分析、控制和预报，好从原始变量所构成的子集合中选择最佳变量，构成最佳变量集合。用主成分分析筛选变量，可以用较少的计算量来选择量，获得选择最佳变量子集合的效果。
#主成分分析法的计算步骤

1. 原始指标数据的标准化采集p维随机变量$x=(X_1,X_2,...,X_p{)}^T,$n个样品$x_i=(x_{i1},x_{i2},...,x_{ip}{)}^T,i=1,2,...,n$
   构造矩阵，对样本矩阵进行标准化变换：
   $$Z_{ij}=\frac{x_{ij}-\overline{x_j}}{s_j},i=1,2,...,n;j=1,2,...,p$$
   其中，$ \bar{x_j}=\frac{\Sigma_{i=1}^n x_{ij}}{n},s_j^{2=\frac{\Sigma_{i=1}}n(x_{ij}-\bar{x_j})^2}{n-1}$得到标准矩阵Z。
2. 对标准化矩阵Z求相关系数矩阵
   $$R=\frac{Z^{T}Z}{n-1}$$
3. 求解样本相关矩阵R的特征向量和特征值（一般软件中是进行奇异值分解）
   相关矩阵R的特征方程$|R-\lambda I_p|=0$得到p个特征根，把p个特征值从大到小排序，求解相应的特征向量。确定需要的主成分个数m，一般使信息的利用率达85%以上。
4. 将标准化后的数据矩阵与m个特征向量相乘，得到m个主成分。
   #R语言实现

假设一组数据：

dat <- data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1),
y=c(2.4,0.7,2.9,2.2,3,2.7,1.6,1.1,1.6,0.9))
dat

   x   y
1  2.5 2.4
2  0.5 0.7
3  2.2 2.9
4  1.9 2.2
5  3.1 3.0
6  2.3 2.7
7  2.0 1.6
8  1.0 1.1
9  1.5 1.6
10 1.1 0.9

数据进行标准化

dat1 <- scale(dat) ##0 mean, 1 std
dat1

               x          y
 [1,]  0.8787452  0.5788568
 [2,] -1.6683424 -1.4294219
 [3,]  0.4966821  1.1695270
 [4,]  0.1146189  0.3425887
 [5,]  1.6428715  1.2876611
 [6,]  0.6240365  0.9332589
 [7,]  0.2419733 -0.3662155
 [8,] -1.0315705 -0.9568857
 [9,] -0.3947986 -0.3662155
[10,] -0.9042161 -1.1931538
attr(,"scaled:center")
   x    y 
1.81 1.91 
attr(,"scaled:scale")
        x         y 
0.7852105 0.8464960 

求解协方差矩阵

covdat <- cov(dat1)
covdat

          x         y
x 1.0000000 0.9259293
y 0.9259293 1.0000000

求解特征值特征向量

eigen_dat <- eigen(covdat)
eigen_dat

$values
[1] 1.92592927 0.07407073

$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

使用奇异值分解时：
$$A=UD{V}^T$$

(a <- svd(covdat))
$d
[1] 1.92592927 0.07407073

$u
           [,1]       [,2]
[1,] -0.7071068 -0.7071068
[2,] -0.7071068  0.7071068

$v
           [,1]       [,2]
[1,] -0.7071068 -0.7071068
[2,] -0.7071068  0.7071068

a$u %*% diag(a$d) %*% t(a$v)

   x         y
x 1.0000000 0.9259293
y 0.9259293 1.0000000

求解主成分

第一个特征值最大，其对应的特征向量为$(0.7071068,0.7071058{)}^T$,现在先算出这两个主成分：

dat1 %*% eigen_dat$vectors

             [,1]        [,2]
 [1,]  1.03068029 -0.21205314
 [2,] -2.19045016  0.16894230
 [3,]  1.17818776  0.47577321
 [4,]  0.32329464  0.16119898
 [5,]  2.07219947 -0.25117173
 [6,]  1.10117414  0.21865330
 [7,] -0.08785251 -0.43005447
 [8,] -1.40605089  0.05281009
 [9,] -0.53811824  0.02021127
[10,] -1.48306451 -0.20430982

用R包自带的函数计算PCA：

pca <- prcomp(dat1)
pca$x

              PC1         PC2
 [1,] -1.03068029  0.21205314
 [2,]  2.19045016 -0.16894230
 [3,] -1.17818776 -0.47577321
 [4,] -0.32329464 -0.16119898
 [5,] -2.07219947  0.25117173
 [6,] -1.10117414 -0.21865330
 [7,]  0.08785251  0.43005447
 [8,]  1.40605089 -0.05281009
 [9,]  0.53811824 -0.02021127
[10,]  1.48306451  0.20430982

pca$rotation

         PC1        PC2
x -0.7071068  0.7071068
y -0.7071068 -0.7071068

可以看出，pca的rotation就是特征向量

把过程整合到一个函数中：

pca <- function(data = data){
  dat <- scale(data)
  covdat <- cov(dat)
  eigendat <- eigen(covdat)
  eigenValue <- eigendat$values
  eigenVector <- eigendat$vectors
  order_value <- order(eigenValue,decreasing = T)
  values <- eigenValue[order_value]
  valueSum <- sum(values)
  cumVar <- cumsum(values)/valueSum * 100
  order_vector <- eigenVector[,order_value]
  principal <- dat %*% order_vector
  return(list(PCA=principal, cumVar=cumVar))
}

pca(data = dat)

$PCA
             [,1]        [,2]
 [1,]  1.03068029 -0.21205314
 [2,] -2.19045016  0.16894230
 [3,]  1.17818776  0.47577321
 [4,]  0.32329464  0.16119898
 [5,]  2.07219947 -0.25117173
 [6,]  1.10117414  0.21865330
 [7,] -0.08785251 -0.43005447
 [8,] -1.40605089  0.05281009
 [9,] -0.53811824  0.02021127
[10,] -1.48306451 -0.20430982

$cumVar
[1]  96.29646 100.00000
 

注意：结果已经按照主成分大小排好序了，第一列为最大主成分，第二列次之，等等。cumVar为累计方差。