stats | 数据降维之主成分分析（PCA）

在使用数据建模时，过多的变量会增加计算复杂性，同时也使结果解释变得困难。主成分分析（Principal Components Analysis，PCA）通过对原变量进行线性组合生成新的变量，可以使得纳入模型中的变量数大大减少，从而在保留大部分信息的前提下达到数据降维的目的。

1 理论基础

记样本的数据矩阵为 ，其中 表示矩阵的第 行、第 列元素； 表示矩阵的总行数，即样本数； 表示矩阵的总列数，即变量数； 表示矩阵的第 行； 表示矩阵的第 列。

1.1 从协方差矩阵求主成分

协方差矩阵

样本矩阵的协方差矩阵 为：

可将 记作 。

特征向量与特征根

记协方差矩阵 的特征根由小到大排列分别为 、 、...、 ，对应的单位特征向量分别为 、 、...、 ，这些特征向量彼此正交。

主成分矩阵

样本矩阵的第 个主成分为：

记 ，则主成分矩阵 。

贡献率

表示第 个主成分的贡献率。

1.2 从相关系数矩阵求主成分

样本矩阵的相关系数矩阵 为：

其中， 为协方差矩阵中的元素。

同样，按照类似于1.1中的步骤，将协方差矩阵 替换成相关系数矩阵 ，也可以得到主成分，称为标准化样本主成分。

2 R中的相关函数

示例数据：

data = USArrests

2.1 cov()/cor()

cov()函数可计算协方差矩阵：

cov(data)
##              Murder   Assault   UrbanPop      Rape
## Murder    18.970465  291.0624   4.386204  22.99141
## Assault  291.062367 6945.1657 312.275102 519.26906
## UrbanPop   4.386204  312.2751 209.518776  55.76808
## Rape      22.991412  519.2691  55.768082  87.72916

cov()函数可计算相关系数矩阵：

cor(data)
##              Murder   Assault   UrbanPop      Rape
## Murder   1.00000000 0.8018733 0.06957262 0.5635788
## Assault  0.80187331 1.0000000 0.25887170 0.6652412
## UrbanPop 0.06957262 0.2588717 1.00000000 0.4113412
## Rape     0.56357883 0.6652412 0.41134124 1.0000000

2.2 eigen()

eigen()函数用于计算矩阵的特征根和特征向量：

eigen(cov(data))
## eigen() decomposition
## $values
## [1] 7011.114851  201.992366   42.112651    6.164246
## 
## $vectors
##             [,1]        [,2]        [,3]        [,4]
## [1,] -0.04170432  0.04482166  0.07989066  0.99492173
## [2,] -0.99522128  0.05876003 -0.06756974 -0.03893830
## [3,] -0.04633575 -0.97685748 -0.20054629  0.05816914
## [4,] -0.07515550 -0.20071807  0.97408059 -0.07232502

eigen(cor(data))
## eigen() decomposition
## $values
## [1] 2.4802416 0.9897652 0.3565632 0.1734301
## 
## $vectors
##            [,1]       [,2]       [,3]        [,4]
## [1,] -0.5358995  0.4181809 -0.3412327  0.64922780
## [2,] -0.5831836  0.1879856 -0.2681484 -0.74340748
## [3,] -0.2781909 -0.8728062 -0.3780158  0.13387773
## [4,] -0.5434321 -0.1673186  0.8177779  0.08902432

• 相关系数矩阵的特征根之和恰好为变量数 。

2.3 prcomp()

prcomp()函数是通过协方差矩阵进行主成分分析的，它的语法结构如下：

prcomp(x, retx = TRUE, center = TRUE,
       scale. = FALSE,
       tol = NULL, rank. = NULL, ...)

• retx：是否输出因子载荷矩阵。

在进行主成分分析时，通常需要对样本数据进行标准化，而prcomp()函数的中心化参数center的默认值是TRUE，而标准化参数scale.的默认值为FALSE，使用时一般需设置scale. = TRUE：

pca_1 <- prcomp(data, scale. = T)

输出结果的数据结构如下：

查看因子载荷矩阵：

pca_1$rotation
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995  0.4181809 -0.3412327  0.64922780
## Assault  -0.5831836  0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158  0.13387773
## Rape     -0.5434321 -0.1673186  0.8177779  0.08902432

查看主成分的标准差：

pca_1$sdev
## [1] 1.5748783 0.9948694 0.5971291 0.4164494

查看各样本的主成分得分（即转换后的变量值）：

head(pca_1$x)
##                   PC1        PC2         PC3          PC4
## Alabama    -0.9756604  1.1220012 -0.43980366  0.154696581
## Alaska     -1.9305379  1.0624269  2.01950027 -0.434175454
## Arizona    -1.7454429 -0.7384595  0.05423025 -0.826264240
## Arkansas    0.1399989  1.1085423  0.11342217 -0.180973554
## California -2.4986128 -1.5274267  0.59254100 -0.338559240
## Colorado   -1.4993407 -0.9776297  1.08400162  0.001450164

绘制碎石图：

screeplot(pca_1, type = "l")

如果并不是样本矩阵的所有变量都参与主成分分析，则可以使用下列语法结构中的formula参数指明参与主成分分析的变量：

prcomp(formula, data = NULL, subset, na.action, ...)

示例：

pca_2 <- prcomp(~ Murder + Assault + UrbanPop, 
                data = data, scale = T)
pca_2$rotation
##                 PC1         PC2        PC3
## Murder   -0.6672955  0.30345520  0.6801703
## Assault  -0.6970818  0.06713997 -0.7138411
## UrbanPop -0.2622854 -0.95047734  0.1667309

2.4 princomp()

princomp()函数是通过相关系数矩阵进行主成分分析的，该函数也比prcomp()函数更加常用，它的语法结构如下：

princomp(x, cor = FALSE, scores = TRUE,
         covmat = NULL,
         subset = rep_len(TRUE, nrow(as.matrix(x))),
         fix_sign = TRUE, ...)

使用princomp()时一般需要设置cor = TRUE：

pca_3 <- princomp(data, cor = T)

输出结果的数据结构如下：

查看因子载荷可以使用loadings()函数：

loadings(pca_3)
## 
## Loadings:
##          Comp.1 Comp.2 Comp.3 Comp.4
## Murder    0.536  0.418  0.341  0.649
## Assault   0.583  0.188  0.268 -0.743
## UrbanPop  0.278 -0.873  0.378  0.134
## Rape      0.543 -0.167 -0.818       
## 
##                Comp.1 Comp.2 Comp.3 Comp.4
## SS loadings      1.00   1.00   1.00   1.00
## Proportion Var   0.25   0.25   0.25   0.25
## Cumulative Var   0.25   0.50   0.75   1.00

绘制碎石图：

screeplot(pca_3, type = "l")

与prcomp()函数一样，它也有如下语法结构：

princomp(formula, data = NULL, subset, na.action, ...)

示例：

pca_4 <- princomp(~ Murder + Assault + UrbanPop, 
                  data = data, cor = T)
loadings(pca_4)
## 
## Loadings:
##          Comp.1 Comp.2 Comp.3
## Murder    0.667  0.303  0.680
## Assault   0.697        -0.714
## UrbanPop  0.262 -0.950  0.167
## 
##                Comp.1 Comp.2 Comp.3
## SS loadings     1.000  1.000  1.000
## Proportion Var  0.333  0.333  0.333
## Cumulative Var  0.333  0.667  1.000