用R語言做主成分分析

前言

主成分分析（PCA)是一種經典的線性維歸約方法，基於高維空間的數據尋找主成分將其轉換至低維空間，並保證低維空間下每個維度具有最大方差且不相關性。主成分分析，就是一種降維的分析方法，其考察多個變量間相關性的一種多元統計方法，研究如何通過少數幾個主成分來揭示多個變量間的內部結構，即從原始變量中導出少數幾個主成分，使它們儘可能多地保留原始變量的信息，且彼此間互不相關。

舉個例子來說明：拿到一個生物信息系的本科生期末考試成績單，裡面有三列，一列是對學科的興趣程度，一列是複習時間，還有一列是考試成績。我們知道要學好相關內容，需要有濃厚的興趣，所以第二項與第一項強相關，第三項和第二項也是強相關。那是不是可以合併第一項和第二項呢？

那麼為什麼要降維呢？

（1）多重共線性-預測變量之間存在一定程度的相關性。

多重共線性會導致解空間的不穩定，從而可能導致結果的不連貫。

（2）高維空間本身具有稀疏性。

（3）過多的變量會妨礙查找規律的建立。

（4）僅在變量層面上分析可能會忽略變量之間的潛在聯繫。

例如幾個預測變量的綁定才可以反映數據某一方面特徵。

主成分計算過程及數據模型

1.  將原始數據標準化，用scale（）函數；

2. 求標準化數據的協方差矩陣，用cov（）函數，或求數據的相關矩陣，用cor（）函數；

3. 求協方差矩陣或相關矩陣的特徵值和單位特徵向量，用eigen（）函數，其中$values是按從大到小對應的特徵值，$vectors是對應的單位特徵向量。

4. 主成分分析，用princomp（x，cor..）函數,x是矩陣，cor確定x是否為相關係數矩陣；

5. 確定主成分個數，用screeplot（）函數，用可視化方法來確定主成分個數，選取一個拐點對應的序號；

6. 解釋主成分，用PCA$loadings顯示主成分載荷矩陣，PCA為主成分分析賦值變量；

7. 確定各樣本的主成分得分，用PCA$scores來確定，並根據樣本各主成分的分值來對樣本進行解釋。

1 函數總結

R中作為主成分分析最主要的函數是princomp()函數

princomp()主成分分析 可以從相關陣或者從協方差陣做主成分分析

summary()提取主成分信息 

loadings()顯示主成分分析或因子分析中載荷的內容

predict()預測主成分的值 

screeplot()畫出主成分的碎石圖 

biplot()畫出數據關於主成分的散點圖和原坐標在主成分下的方向

1 案例演示

現有30名中學生身高、體重、胸圍、坐高數據，對身體的四項指標數據做主成分分析。

(1)載入原始數據

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))

（2）做主成分分析並顯示分析結果

test.pr<-princomp(test,cor=TRUE) 

#cor是邏輯變量，當cor=TRUE時表示用樣本的相關矩陣R做主成分分析；當cor=FALSE時表示用樣本的協方差陣S做主成分分析

summary(test.pr,loadings=TRUE) 

#loading是邏輯變量，當loading=TRUE時表示顯示loading 的內容

#loadings的輸出結果為載荷，是主成分對應於原始變量的係數即Q矩陣

（3）分析結果含義

Standard deviation 標準差  其平方為方差=特徵值

Proportion of Variance  方差貢獻率

Cumulative Proportion  方差累計貢獻率

由結果顯示 前兩個主成分的累計貢獻率已經達到96% 可以捨去另外兩個主成分，達到降維的目的

因此可以得到函數表達式：

Z1=-0.497X1-0.515X2-0.481X3-0.507X4

Z1=0.543X1-0.210X2-0.725X3-0.368X4

（4）畫主成分的碎石圖並預測

screeplot(test.pr,type="lines")

由圖可以看出，第二個主成分之後圖線變化趨於平穩，因此可以選擇前兩個主成分做分析

p<-predict(test.pr)

# predict函數根據主成分生成最後的新變量(即將原始數據矩陣X與主成分係數矩陣A結合,生成最後的主成分新變量,用這個生成的新變量代替原始變量)

通訊郵箱：[e-mail]

投稿郵箱：[e-mail]

往期「精彩內容」，點擊回顧

DNA測序歷史 | CircRNA資料庫 | Epigenie表觀綜合 | 癌症定位

BWA介紹  |  源碼安裝R包 | CancerLocator  | lme4 | 450K分析

乳腺癌異質性 |  BS-Seq  | 隱馬模型 | Circos安裝 |  Circos畫圖

KEGG標記基因 |  GDSC  |  Meta分析  |   R線性回歸和相關矩陣

精彩會議及課程，點擊回顧

計算表觀遺傳學大數據前沿學術論壇會議記實

哈爾濱醫科大學2017年全國生物信息學暑期學校

2017龍星課程系列（一）

2017龍星課程系列（二）
2017龍星課程系列（三）

2017龍星課程系列（四）

2017龍星課程系列（五）

CEPI感謝您的支持！

（IOS系統用戶專用通道）

    讚賞