主成分分析

PCA原理 - 奇异值分解

假设有一组包含\(m\)个变量的数据，含有\(n\)个样本，构成了一个\(m \times n\)的矩阵\(A_0\)。用图形表示，\(A_0\)是指在\(R^m\)空间内的\(n\)个点，当\(A\)的每一行都减去该行的均值后（数据中心化），\(n\)个点通常聚集在一条线或者平面周围（或其他\(R^m\)的低维子空间)。

以含两个变量的\(n\)个点为例说明PCA的步骤，\(A\) 维度为\(2 \times n\)。\(A\)的每一行减去该行的均值，使\(A\)中心化，每一行的均值为0。如图所示，这\(n\)个点分布在一条直线周围。由图中也可以直观的看到，主成分方向还表示原始数据在该方向上投影方差最大（分散），与之垂直的投影距离最小。

\(AA^T\)构成一个\(2 \times 2\)的矩阵，其中\((AA^T)_{11}\)和\((AA^T)_{22}\)是行向量的方差，而\((AA^T)_{12}\)和\((AA^T)_{21}\)即是两行的协方差。因此，协方差矩阵可表示为

\[ S = \frac {AA^T} {n-1} \]

中心化的原因：保证第一个主分量是方差最大的方向，如果不中心化，那么第一主成分的方向受均值的影响而有一定的倾向性，如下图。

求协方差矩阵的特征值\(\sigma_1\)和\(\sigma_1\)。根据矩阵奇异值分解公式，中心化后的\(A\)可表示为

\[ A = \sqrt{\sigma_1}u_1v_1^T + \sqrt{\sigma_2}u_2v_2^T \]

假设中心化后\(A\)为（各行和为0）

\[ A = \left[ \begin{matrix} 3 & -4 & 7 & 1 & -4 & -3 \\ 7 & -6 & 8 & -1 & -1 & -7 \end{matrix} \right] \qquad S = \frac{AA^T}{5} = \left[ \begin{matrix} 20 & 25 \\ 25 & 40 \end{matrix} \right] \]

\(S\)特征值为57和3，那么\(A\)可表示为

\[ A = \sqrt{57}u_1v_1^T + \sqrt{3}u_2v_2^T \]

\(u_1\)对应于散点图变化最大的方向（直线），\(u_2\)垂直于\(u_1\)，表征垂直于直线方向上的小幅摆动。

所以，对于多维变量，把原始数据转换到向量\(u\)构成的空间内，特征值越大表征在该向量的方向上数据变异越大。通常选择变异量较大的前几个变量用于表征数据。

PCA数学推导

以最大投影方差法推导，原始数据的中心点为：

\[ \bar{x} = \frac{\sum_{n=1}^{N}}{N} \] 原始数据向投影向量\(u\)投影之后的方差为：

\[ \frac{1}{N} \sum_{n-1}^{N}(u^Tx_n - u^T\bar{x})^2 = u^TSu \] 根据梯度优化原则（拉格朗日乘子）:

\[ u^TSu + \lambda(1-u^Tu) = 0 \] \[ Su = \lambda u \]

因此，对于多维数据，协方差矩阵\(S\)的特征值即为方差最大方向，从公式中也可看出PCA分析需要预先对数据中心化。

library(dplyr)
# 首先中心化，因为各变量量纲不同scale = TRUE,
# 主成分是原始变量的线性组合，$rotation表征了线性组合系数，x
# 表示新的坐标值
mtcars_pca <- prcomp(select(mtcars, -c("vs", "am")),
  center = TRUE, scale = TRUE)

summary(mtcars_pca)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6
## Standard deviation     2.3782 1.4429 0.71008 0.51481 0.42797 0.35184
## Proportion of Variance 0.6284 0.2313 0.05602 0.02945 0.02035 0.01375
## Cumulative Proportion  0.6284 0.8598 0.91581 0.94525 0.96560 0.97936
##                            PC7    PC8     PC9
## Standard deviation     0.32413 0.2419 0.14896
## Proportion of Variance 0.01167 0.0065 0.00247
## Cumulative Proportion  0.99103 0.9975 1.00000

library(ggbiplot)
ggbiplot(mtcars_pca, labels=rownames(mtcars))

mtcars_country <- c(rep("Japan", 3), rep("US",4), rep("Europe", 7),rep("US",3), "Europe", rep("Japan", 3), rep("US",4), rep("Europe", 3), "US", rep("Europe", 3))

ggbiplot(mtcars_pca,ellipse=TRUE,  labels=rownames(mtcars), groups=mtcars_country)

ggbiplot(mtcars_pca,ellipse=TRUE,choices=c(3,4),   labels=rownames(mtcars), groups=mtcars_country)