为什么要做主成分分析

不管三七二十一就直接套用统计方法都是耍流氓，做主成分分析并不是拍脑袋决定的。 在这个例子里面，我们拿到了这43个法官的12个信息，就可以通过这12个指标来对法官进行分类，但也许实际情况下收集其他法官的12个信息比较麻烦，所以我们希望只收集三五个信息即可，然后也想达到比较好的分类效果。或者至少也得剔除几个指标吧，这个时候主成分分析就能派上用场啦。这12个变量能得到12个主成分，如果前两个主成分可以揭示85%以上的变异度，也就是说我们可以用两个主成分来代替那12个指标。

在生物信息学领域，比如我们测了1000个病人的2万个基因的表达矩阵，同时也有他们的健康状态信息。那么我们想仔细研究这些数据，想得到基因表达与健康状态的某种关系。这样我就可以对其余几十亿的人检测基因表达来预测其健康状态。如果我们进行了主成分分析，就可以选择解释度比较高的主成分对应的基因，可能就几十上百个而已，大幅度的降低广泛的基因检测成本。

• 数据标准化(中心化)

  dat_scale=scale(USJudgeRatings,scale=F)
  options(digits=4, scipen=4)
  kable(head(dat_scale)) #scale()是对数据中心化的函数，当参数scale=F时，表示将数据按列减去平均值，scale=T表示按列进行标准化，公式为（x-mean(x)）/sd(x)，本例采用中心化。

• 求相关系数矩阵

  dat_cor=cor(dat_scale)#相关系数，协方差是cov
  options(digits=4, scipen=4)
  kable(dat_cor)

• r语言options

digits ：指定在输出数字时要显示的小数位数。
scipen ：控制科学计数法中超过多少位数时切换为科学计数法。
warn ：指定在发生警告时是否输出警告信息。
stringsAsFactors ：控制是否将字符串类型的变量自动转换为因子类型。
max.print ：指定在输出时最多显示多少行

• 计算特征值和特征向量

  dat_eigen=eigen(dat_cor)
  dat_var=dat_eigen$values ## 相关系数矩阵的特征值
  options(digits=4, scipen=4)
  dat_var
  ##  [1] 10.133504  1.104147  0.332902  0.253847  0.084453  0.037286  0.019683
  ##  [8]  0.015415  0.007833  0.005612  0.003258  0.002060
  pca_var=dat_var/sum(dat_var)
  pca_var
  ##  [1] 0.8444586 0.0920122 0.0277418 0.0211539 0.0070377 0.0031072 0.0016402
  ##  [8] 0.0012846 0.0006528 0.0004676 0.0002715 0.0001717
  pca_cvar=cumsum(dat_var)/sum(dat_var)
  pca_cvar
  ##  [1] 0.8445 0.9365 0.9642 0.9854 0.9924 0.9955 0.9972 0.9984 0.9991 0.9996
  ## [11] 0.9998 1.0000

• 崖低碎石图和累积贡献图

  library(ggplot2)
  p=ggplot(,aes(x=1:12,y=pca_var))
  p1=ggplot(,aes(x=1:12,y=pca_cvar))
  p+geom_point(pch=2,lwd=3,col=2)+geom_line(col=2,lwd=1.2)

p1+geom_point(pch=2,lwd=3,col=2)+geom_line(col=2,lwd=1.2)

• 主成分载荷

  pca_vect= dat_eigen$vector  ## 相关系数矩阵的特征向量
  loadings=sweep(pca_vect,2,sqrt(pca_var),"*")
  rownames(loadings)=colnames(USJudgeRatings)
  options(digits=4, scipen=4)
  kable(loadings[,1:2])