主成分分析（PCA）及案例分析

主成分分析（parincipal component analysis,PCA）

#对 USA ests 数据集进行 PCA， PCA 包肯在基础软件包中。

states=row.names(USArrests)  #数据集包含50个州

states #显示50个州的名字

 

names(USArrests)  #数据集的列包含4个变量

 

apply(USArrests, 2, mean) #计算数据集每列的均值，1表示对行操作，2表示对列操作

 

结果分析：数据显示强奸案平均发生次数是谋杀案的3倍，袭击案平均发生次数是强奸案的 8倍以上，变量的均值差异很大。

apply(USArrests, 2, var) #apply ()函数计算4个变量的方差。

 

结果分析：变量的方差之间存在着较大的差异，。如果不对变量标准化就进行 PCA ，就会导致大多数主成分会出 Assault 一个变量所决定，因为 Assault变量的均值和方差明显是最大的。因此，在进行 PCA 之前对变最进行标准化处理是非常必要的。

pr.out=prcomp(USArrests, scale=TRUE) #prcomp ()函数进行主成分分析，prcomp ()函数默认对变量进行中心化处理，选项 scale = TURE 可以对变最进行标准化处理。

names(pr.out)

 

结果分析：center和 scale 表示在实施 PCA 之前进行标准化以后变量的均值和标准差，rotation 矩阵提供了主成分载荷信息。

pr.out$center

 

pr.out$scale

 

pr.out$rotation #pr. Out$rotation的列包含对应的主成分载荷向量

 

结果分析：可以看到有4个不同的主成分，因为在一个有n个观测和p个变量的数据集中一般有min(n-1,p)个信息量在较大的主成分。第二主成分向量PC2在 UrbanPop 上有较大权重，在其他3个变量上权重较小，因此，这个主成分大致刻画了每个州的城市化水平。

dim(pr.out$x) #主成分得分向量长度n=50，主成分载荷向量的长度为 =4

 

biplot(pr.out, scale=0) #绘制前两个主成分的双标图，参数scale =0 表示载荷箭头所指的变量经过了标准化

 

结果分析：第一载荷向量在 Assault 、Murder 、Rape这3个变量上的权重大致相等，而在 UrbanPop上的权重则相对较小，因此这个主成分大致解释了严重罪行的总体犯罪率。总之，与犯罪相关的变盘 (Assault、Murder和Rape) 之间的位置比较接近， Urbanpop变量与这3者相对较远。这表明与犯罪相关的变量彼此相关性较强，即谋杀率高的州抢劫和强奸的比率也会比较高，变量UrbanPop 与其他三者相关性较弱。

#主成分有一个性质是:在符号可变的意义下唯一，所以可以通过一些小的改变重新绘制

pr.out$rotation=-pr.out$rotation

pr.out$x=-pr.out$x

biplot(pr.out, scale=0)

 

pr.out$sdev #输出每个主成分的标准差

 

pr.var=pr.out$sdev^2 #主成分解释的方差

pr.var

 

pve=pr.var/sum(pr.var) #要计算每个主成分的方差解释比例，用每个主成分解释的方差除以 4个主成分解释的总方差。

pve

 

结果分析：可以看到第一主成分解释了数据中 62% 的方差，第二主成分解释了数据中 24.7% 的方差。

plot(pve,xlab="Principal Component",ylab="Proportion of Variance Explained", ylim=c(0,1),type='b') #以绘制每个主成分的 PVE 图

 

plot(cumsum(pve), xlab="Principal Component", ylab="Cumulative Proportion of Variance Explained", ylim=c(0,1),type='b') #绘制每个主成分累积 PVE 的图，cumsum ()函数用于计算数值向量中的元素的累计和

 

#计算累计和的例子

a=c(1,2,8,-3)

cumsum(a) 

 

案例1：NCI60数据集的应用

#NCI60数据集由64个细胞系的6830个基因表达数据构成，每个细胞系都有一个标签变量记录了其癌细胞的类型

library(ISLR)

nci.labs=NCI60$labs

nci.data=NCI60$data

#在进行 PCA 和聚类分析时，由于是无指导约学习，不使用癌细胞的类型特征，在 PCA和聚类分析之后，可以用这些特征检验无指导学习的结果与癌细胞类型的匹配度。

dim(nci.data)

 

nci.labs[1:4] #查看癌细胞系的类型

 

table(nci.labs) #列出癌细胞的类型

 

pr.out=prcomp(nci.data, scale=TRUE) #对变量进行标准化

#编写一个可以给数值向量每个元素分配不同颜色的简单函数，这个函数可以根据每个观测对应的癌症类型给64个细胞系分配不同的颜色。

Cols=function(vec){

  cols=rainbow(length(unique(vec)))

  return(cols[as.numeric(as.factor(vec))])

}

par(mfrow=c(1,2))

plot(pr.out$x[,1:2], col=Cols(nci.labs), pch=19,xlab="Z1",ylab="Z2") #绘制主成分的得分向量图

 

plot(pr.out$x[,c(1,3)], col=Cols(nci.labs), pch=19,xlab="Z1",ylab="Z3") ##绘制主成分的得分向量图

 

结果分析：。从图上来看，对应于同类癌症的细胞在前几个主成分得分向量上的值确实更接近，这表明间一类癌症的细胞系往往有非常相似的基因表达水平。

summary(pr.out) #用summary ()函数得到前几个主成分的方差解释比例 (PVE) 的汇总表

 

plot(pr.out) #数绘制出前几个主成分解释的方差的图

 

结果分析：，在柱形图中，每个柱子的高度是 pr. out $ sdev 相应元索的平方。

pve=100*pr.out$sdev^2/sum(pr.out$sdev^2)

par(mfrow=c(1,2))

plot(pve,  type="o", ylab="PVE", xlab="Principal Component", col="blue")

 

plot(cumsum(pve), type="o", ylab="Cumulative PVE", xlab="Principal Component", col="brown3")

 

结果分析：前7个主成分一共解释了数据大约 40% 的方差， 但这个比例还不够大。然而，通过碎石图观察发现前7个卖成分中每一个都解释了大量方差，而之后的主成分对方差的解释作用明显下降，即大约在碎石图中的第7个主成分的位置量有一个肘。这表示没有必要考虑7个以上的主成分(尽管对7个主成分分析也已经很有挑战了)。

 

#可以画碎石图，帮助判断主成分个数

screeplot(pca,type="line",main="碎石图")