主成份分析（Principal component analysis, PCA） - ricket-sjtu/bi028 GitHub Wiki

主成份分析的问题

作为一种数据降维的手段，需要用到几个主成份？
主成份得分的意义是什么？
主成份分析中，什么时候该用variance-covariance矩阵？什么时候用correlation矩阵？

PCA工作流程

计算样本方差-协方差矩阵$\textbf{S} = \frac{1}{n-1} \sum_{i=1}^{n}(\mathbf{X}_i-\bar{\textbf{x}})(\mathbf{X}_i-\bar{\textbf{x}})'$
计算$\mathbf{S}$的特征值$\hat{\lambda}_1, \hat{\lambda}_2, \ldots, \hat{\lambda}_p$，以及对应的特征向量$\hat{\mathbf{e}}_1, \hat{\mathbf{e}}_2, \ldots, \hat{\mathbf{e}}_p$
计算主成份：$\hat{\mathbf{Y}}i = \hat{\mathbf{e}}{i1}\mathbf{X}1 + \hat{\mathbf{e}}{i2}\mathbf{X}2 + \ldots + \hat{\mathbf{e}}{ip}\mathbf{X}_p$

主成份的选择

一般来说，我们需要前$k$个主成份：

$k$必须足够小，以便于理解；
$k$必须足够大，能够解释足够比例的变化（variation）： $$ \frac{\hat{\lambda}_1 + \hat{\lambda}_2 + \dots + \hat{\lambda}_k}{\hat{\lambda}_1 + \hat{\lambda}_2 + \dots + \hat{\lambda}_p} \cong 1 $$

如何选择用方差-协方差矩阵还是相关矩阵

如果使用原始数据的方差-协方差矩阵，前几个主成份中，原始数据中方差较大的几个变量将占主导；
如果原始数据的不同变量具有相同的单位，则采用原始数据的vcov矩阵；
如果原始数据的不同变量单位不同，必须对每个变量进行标准化处理后再计算vcov矩阵，这时候计算得到的vcov矩阵相当于correlation矩阵；
因此，选择用原始的vcov矩阵还是correlation矩阵取决于是否给不同变量相同的权重（correlation矩阵）还是不同的权重（raw vcov矩阵）

PCA Demo

R中作为主成分分析最主要的函数是princomp()函数
- princomp()可以从相关阵或者从协方差阵做主成分分析
- summary()提取主成分信息
- loadings()显示主成分分析或因子分析中载荷的内容
- predict()预测主成分的值
- screeplot()画出主成分的碎石图
- biplot()画出数据关于主成分的散点图和原坐标在主成分下的方向
案例

test<-data.frame(
  X1=c(148, 139, 160, 149, 159, 142, 153, 150, 151, 139,
           140, 161, 158, 140, 137, 152, 149, 145, 160, 156,
           151, 147, 157, 147, 157, 151, 144, 141, 139, 148),
  X2=c(41, 34, 49, 36, 45, 31, 43, 43, 42, 31,
           29, 47, 49, 33, 31, 35, 47, 35, 47, 44,
           42, 38, 39, 30, 48, 36, 36, 30, 32, 38),
  X3=c(72, 71, 77, 67, 80, 66, 76, 77, 77, 68,
          64, 78, 78, 67, 66, 73, 82, 70, 74, 78,
          73, 73, 68, 65, 80, 74, 68, 67, 68, 70),
  X4=c(78, 76, 86, 79, 86, 76, 83, 79, 80, 74,
           74, 84, 83, 77, 73, 79, 79, 77, 87, 85,
           82, 78, 80, 75, 88, 80, 76, 76, 73, 78)
  )
#cor是逻辑变量 当cor=TRUE表示用样本的相关矩阵R做主成分分析
#当cor=FALSE表示用样本的协方差阵S做主成分分析
test.pr<-princomp(test,cor=TRUE)
#loading是逻辑变量 当loading=TRUE时表示显示loading 的内容
#loadings的输出结果为载荷 是主成分对应于原始变量的系数即Q矩阵
summary(test.pr,loadings=TRUE)

结果显示

Importance of components:
                          Comp.1     Comp.2     Comp.3     Comp.4
Standard deviation     1.8817805 0.55980636 0.28179594 0.25711844
Proportion of Variance 0.8852745 0.07834579 0.01985224 0.01652747
Cumulative Proportion  0.8852745 0.96362029 0.98347253 1.00000000

Loadings:
   Comp.1 Comp.2 Comp.3 Comp.4
X1 -0.497  0.543 -0.450  0.506
X2 -0.515 -0.210 -0.462 -0.691
X3 -0.481 -0.725  0.175  0.461
X4 -0.507  0.368  0.744 -0.232

结果的意义

Importance of components：主成份的重要性度量
Standard deviation：标准偏差，其平方是方差，也是对应的特征值；
Proportion of Variance：方差占总方差的比例；
Cumulative Proportion：累积占总方差的比例，前两个主成份的累积贡献率已经达到96.4%，可以舍去其他两个主成份，也就是说$k=2$；
Loadings：载荷矩阵，也就是对应于每个主成份的变量线性组合的系数
- Comp1 = -0.497*X1 - 0.515*X2 - 0.481*X3 - 0.507*X4

PCA结果的绘图表示

# 碎石图screeplot
screeplot(test.pr) # 贡献率默认用barplot表示
screeplot(test.pr,type="lines")  # 贡献率的线型表示
# 散点图，只绘出前两个主成份，以及变量对应的向量
biplot(test.pr)

第二种PCA

分析也可从相关系数矩阵出发：

# 1. 下三角相关系数矩阵，先用向量表示
x<-c(1.00,
        0.79, 1.00,
        0.36, 0.31, 1.00,
        0.96, 0.74, 0.38, 1.00,
        0.89, 0.58, 0.31, 0.90, 1.00,
        0.79, 0.58, 0.30, 0.78, 0.79, 1.00,
        0.76, 0.55, 0.35, 0.75, 0.74, 0.73, 1.00,
        0.26, 0.19, 0.58, 0.25, 0.25, 0.18, 0.24, 1.00,
        0.21, 0.07, 0.28, 0.20, 0.18, 0.18, 0.29,-0.04, 1.00,
        0.26, 0.16, 0.33, 0.22, 0.23, 0.23, 0.25, 0.49,-0.34, 1.00,
        0.07, 0.21, 0.38, 0.08,-0.02, 0.00, 0.10, 0.44,-0.16, 0.23, 1.00,
        0.52, 0.41, 0.35, 0.53, 0.48, 0.38, 0.44, 0.30,-0.05, 0.50,0.24, 1.00,
        0.77, 0.47, 0.41, 0.79, 0.79, 0.69, 0.67, 0.32, 0.23, 0.31,0.10, 0.62, 1.00,
        0.25, 0.17, 0.64, 0.27, 0.27, 0.14, 0.16, 0.51, 0.21, 0.15,0.31, 0.17, 0.26, 1.00,
        0.51, 0.35, 0.58, 0.57, 0.51, 0.26, 0.38, 0.51, 0.15, 0.29,0.28, 0.41, 0.50, 0.63, 1.00,
        0.21, 0.16, 0.51, 0.26, 0.23, 0.00, 0.12, 0.38, 0.18, 0.14,0.31, 0.18, 0.24, 0.50, 0.65, 1.00)
# 2.定义变量名
names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8",
                  "X9","X10", "X11", "X12", "X13", "X14", "X15", "X16")
# 3.转化为相关系数矩阵
R<-matrix(0,nrow=16,ncol=16,dimnames=list(names,names))
for(i in 1:16){
  for(j in 1:i){
    R[i,j]<-x[(i-1)*i/2+j];R[j,i]<-R[i,j]
  }
}
# 4.主成分分析
x.pr<-princomp(covmat=R)  #covmat协方差阵
load<-loadings(x.pr)  #载荷

# 5.散点图
plot(load[,1:2])
text(load[,1],load[,2],adj=c(-0.4,0.3))

可以根据loadings对变量进行聚类，当然也可以直接在散点图上对变量进行聚类。