主成分分析（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="碎石圖")