• 1 背景

• 2 拆解主成分分析步驟

• 2.1 測(cè)試數(shù)據(jù)

• 2.2 為什么要做主成分分析

• 2.3 step1:數(shù)據(jù)標(biāo)準(zhǔn)化(中心化)

• 2.4 step2:求相關(guān)系數(shù)矩陣

• 2.5 step3:計(jì)算特征值和特征向量

• 2.6 step4:崖低碎石圖和累積貢獻(xiàn)圖

• 2.7 step5:主成分載荷

• 2.8 step6:主成分得分計(jì)算和圖示

• 3 實(shí)戰(zhàn)一

• 4 實(shí)戰(zhàn)二

• 5 進(jìn)階的主成分分析-psych包

• 6 推薦一個(gè)R包factoextra

• 7.主成分分析的生物信息學(xué)應(yīng)用

• 8. 主成分分析的其它可視化方法

• 9.其它學(xué)習(xí)資料
  

1 背景

主成分分析法是數(shù)據(jù)挖掘中常用的一種降維算法,是Pearson在1901年提出的,再后來由hotelling在1933年加以發(fā)展提出的一種多變量的統(tǒng)計(jì)方法，其最主要的用途在于“降維”，通過析取主成分顯出的最大的個(gè)別差異,也可以用來削減回歸分析和聚類分析中變量的數(shù)目，與因子分析類似。

所謂降維，就是把具有相關(guān)性的變量數(shù)目減少，用較少的變量來取代原先變量。如果原始變量互相正交，即沒有相關(guān)性，則主成分分析沒有效果。

在生物信息學(xué)的實(shí)際應(yīng)用情況中，通常是得到了成百上千個(gè)基因的信息，這些基因相互之間會(huì)有影響，通過主成分分析后，得到有限的幾個(gè)主成分就可以代表它們的基因了。也就是所謂的降維。

R語言有非常多的途徑做主成分分析，比如自帶的princomp()和psych包的principal()函數(shù)，還有g(shù)models包的fast.prcomp函數(shù)。

2 拆解主成分分析步驟

實(shí)際應(yīng)用時(shí)我們通常會(huì)選擇主成分分析函數(shù)，直接把input數(shù)據(jù)一步分析到位，只需要看懂輸出結(jié)果即可。但是為了加深理解，這里一步步拆解主成分分析步驟，講解原理。

2.1 測(cè)試數(shù)據(jù)

數(shù)據(jù)集USJudgeRatings包含了律師對(duì)美國(guó)高等法院法官的評(píng)分。數(shù)據(jù)框包含43個(gè)樣本，12個(gè)變量!

下面簡(jiǎn)單看一看這12個(gè)變量是什么，以及它們的相關(guān)性。

這12個(gè)變量的介紹如下：

這些是專業(yè)領(lǐng)域的用詞，大家可以先不用糾結(jié)具體細(xì)節(jié)。

2.2 為什么要做主成分分析

不管三七二十一就直接套用統(tǒng)計(jì)方法都是耍流氓，做主成分分析并不是拍腦袋決定的。  在這個(gè)例子里面，我們拿到了這43個(gè)法官的12個(gè)信息，就可以通過這12個(gè)指標(biāo)來對(duì)法官進(jìn)行分類，但也許實(shí)際情況下收集其他法官的12個(gè)信息比較麻煩，所以我們希望只收集三五個(gè)信息即可，然后也想達(dá)到比較好的分類效果?；蛘咧辽僖驳锰蕹龓讉€(gè)指標(biāo)吧，這個(gè)時(shí)候主成分分析就能派上用場(chǎng)啦。這12個(gè)變量能得到12個(gè)主成分，如果前兩個(gè)主成分可以揭示85%以上的變異度，也就是說我們可以用兩個(gè)主成分來代替那12個(gè)指標(biāo)。

在生物信息學(xué)領(lǐng)域，比如我們測(cè)了1000個(gè)病人的2萬個(gè)基因的表達(dá)矩陣，同時(shí)也有他們的健康狀態(tài)信息。那么我們想仔細(xì)研究這些數(shù)據(jù)，想得到基因表達(dá)與健康狀態(tài)的某種關(guān)系。這樣我就可以對(duì)其余幾十億的人檢測(cè)基因表達(dá)來預(yù)測(cè)其健康狀態(tài)。如果我們進(jìn)行了主成分分析，就可以選擇解釋度比較高的主成分對(duì)應(yīng)的基因，可能就幾十上百個(gè)而已，大幅度的降低廣泛的基因檢測(cè)成本。

2.3 step1:數(shù)據(jù)標(biāo)準(zhǔn)化(中心化)

scale()是對(duì)數(shù)據(jù)中心化的函數(shù)，當(dāng)參數(shù)scale=F時(shí)，表示將數(shù)據(jù)按列減去平均值，scale=T表示按列進(jìn)行標(biāo)準(zhǔn)化，公式為（x-mean(x)）/sd(x)，本例采用中心化。

2.4 step2:求相關(guān)系數(shù)矩陣

從相關(guān)系數(shù)看，只有 CONT 變量跟其它變量沒有相關(guān)性。

當(dāng)然，這樣的矩陣不易看清楚規(guī)律，很明顯，畫個(gè)熱圖就一眼看出來了。

2.5 step3:計(jì)算特征值和特征向量

利用eigen函數(shù)計(jì)算相關(guān)系數(shù)矩陣的特征值和特征向量。這個(gè)是主成分分析法的精髓。

可以看出，PC1(84.4%)和PC2(9.2%)共可以解釋這12個(gè)變量的93.6的程度，除了CONT外的其他的11個(gè)變量與PC1都有較好的相關(guān)性，所以PC1與這11個(gè)變量基本斜交，而CONT不能被PC1表示，所以基本與PC1正交垂直，而PC2與CONT基本平行，表示其基本可以表示CONT。

2.6 step4:崖低碎石圖和累積貢獻(xiàn)圖

崖低碎石圖（scree plot）即貢獻(xiàn)率圖，是希望圖形一開始很陡峭，如懸崖一般，而剩下的數(shù)值都很小，如崖底的碎石一樣。

崖低碎石圖和累積貢獻(xiàn)率圖是對(duì)主成分貢獻(xiàn)率和累積貢獻(xiàn)率的一種直觀表示，用以作為選擇主成分個(gè)數(shù)的參考。本例中第一個(gè)主成分解釋總變異的84.4%,可以只選擇第一個(gè)主成分，但第二主成分也不小，因此選擇前兩個(gè)主成分。

主成分的個(gè)數(shù)選擇沒有一定之規(guī)，需按實(shí)際情況具體分析，一般要求累積貢獻(xiàn)率大于85%或特征值大于1.

但是在實(shí)際的生物信息學(xué)應(yīng)用中，通常達(dá)不到這個(gè)要求。

2.7 step5:主成分載荷

主成分載荷表示各個(gè)主成分與原始變量的相關(guān)系數(shù)。

結(jié)果表明，CONT指標(biāo)跟其它指標(biāo)表現(xiàn)完全不一樣，第一個(gè)主成分很明顯跟除了CONT之外的所有其它指標(biāo)負(fù)相關(guān)，而第二個(gè)主成分則主要取決于CONT指標(biāo)。

2.8 step6:主成分得分計(jì)算和圖示

將中心化的變量dat_scale乘以特征向量矩陣即得到每個(gè)觀測(cè)值的得分。

將兩個(gè)主成分以散點(diǎn)圖形式展示,看看這些樣本被這兩個(gè)主成分是如何分開的

對(duì)于主成分分析，不同數(shù)據(jù)會(huì)有不同的分析方法，應(yīng)具體情況具體分析。

總結(jié)一下PCA的算法步驟：

設(shè)有m條n維數(shù)據(jù)。

• 1）將原始數(shù)據(jù)按列組成n行m列矩陣X

• 2）將X的每一行（代表一個(gè)屬性字段）進(jìn)行零均值化，即減去這一行的均值

• 3）求出協(xié)方差矩陣

• 4）求出協(xié)方差矩陣的特征值及對(duì)應(yīng)的特征向量

• 5）將特征向量按對(duì)應(yīng)特征值大小從上到下按行排列成矩陣，取前k行組成矩陣P

• 6）Y=PX即為降維到k維后的數(shù)據(jù)

PCA本質(zhì)上是將方差最大的方向作為主要特征，并且在各個(gè)正交方向上將數(shù)據(jù)“離相關(guān)”，也就是讓它們?cè)诓煌环较蛏蠜]有相關(guān)性。

PCA也存在一些限制，例如它可以很好的解除線性相關(guān)，但是對(duì)于高階相關(guān)性就沒有辦法了，對(duì)于存在高階相關(guān)性的數(shù)據(jù)，可以考慮Kernel  PCA，通過Kernel函數(shù)將非線性相關(guān)轉(zhuǎn)為線性相關(guān)，關(guān)于這點(diǎn)就不展開討論了。另外，PCA假設(shè)數(shù)據(jù)各主特征是分布在正交方向上，如果在非正交方向上存在幾個(gè)方差較大的方向，PCA的效果就大打折扣了。

最后需要說明的是，PCA是一種無參數(shù)技術(shù)，也就是說面對(duì)同樣的數(shù)據(jù)，如果不考慮清洗，誰來做結(jié)果都一樣，沒有主觀參數(shù)的介入，所以PCA便于通用實(shí)現(xiàn)，但是本身無法個(gè)性化的優(yōu)化。

3 實(shí)戰(zhàn)一

比如你要做一項(xiàng)分析人的糖尿病的因素有哪些，這時(shí)你設(shè)計(jì)了10個(gè)你覺得都很重要的指標(biāo)，然而這10個(gè)指標(biāo)對(duì)于你的分析確實(shí)太過繁雜，這時(shí)你就可以采用主成分分析的方法進(jìn)行降維。10個(gè)指標(biāo)之間會(huì)有這樣那樣的聯(lián)系，相互之間會(huì)有影響，通過主成分分析后，得到三五個(gè)主成分指標(biāo)。此時(shí),這幾個(gè)主成分指標(biāo)既涵蓋了你10個(gè)指標(biāo)中的絕大部分信息，這讓你的分析得到了簡(jiǎn)化（從10維降到3、5維）。

下面是442個(gè)糖尿病人相關(guān)的數(shù)據(jù)集，具體如下：

• x a matrix with 10 columns (自變量)

• y a numeric vector (因變量)

其中x矩陣含有10個(gè)變量，分別是：“age” “sex” “bmi” “map” “tc” “l(fā)dl” “hdl” “tch” “l(fā)tg” “glu” 它們都在一定程度上或多或少的會(huì)影響個(gè)體糖尿病狀態(tài)。

數(shù)據(jù)的詳細(xì)介紹見 Efron, Hastie, Johnstone and Tibshirani (2003) “Least Angle Regression” (with discussion) Annals of Statistics;

一步法主成分分析

上面我們把整個(gè)主成分分析步驟拆解開來講解具體原理，但是實(shí)際數(shù)據(jù)處理過程中我們通常是用現(xiàn)成的函數(shù)一步法完成主成分分析，而且這個(gè)是非常高頻的分析，所以R里面自帶了一個(gè)函數(shù)princomp()來完成主成分分析，如下：

cor是邏輯變量,當(dāng)cor=TRUE表示用樣本的相關(guān)矩陣R做主成分分析,當(dāng)cor=FALSE表示用樣本的協(xié)方差陣S做主成分分。

可以看到前三個(gè)主成份的信息量也只有67.2%，達(dá)不到我們前面說到85%，所以很難說可以用這3個(gè)主成分去代替這10個(gè)生理指標(biāo)來量化病人的狀態(tài)。

值得一提的是，如果你看懂了前面的主成分分析的拆解步驟，就應(yīng)該明白有多少個(gè)變量就有多少個(gè)主成分，只是并不是所有的主成分都有意義，理想狀態(tài)下我們希望有限的幾個(gè)主成分就可以代替數(shù)量繁多的變量，尤其是生物信息學(xué)里面的基因表達(dá)矩陣，兩三萬個(gè)基因的表達(dá)情況我們希望幾十個(gè)基因就可以替代它們，因?yàn)槟切┗蛑g是相互關(guān)聯(lián)的。

碎石圖

也可以畫出主成分的碎石圖，來決定選幾個(gè)主成分。

由碎石圖可以看出，第5個(gè)主成分之后，圖線變化趨于平穩(wěn)，因此可以選擇前5個(gè)主成分做分析。

樣本分布的散點(diǎn)圖

根據(jù)前兩個(gè)主成分畫出樣本分布的散點(diǎn)圖。

上面這個(gè)圖似乎意義不大，因?yàn)榇蟛糠智闆r下都是需要結(jié)合樣本的分組信息來看看這些主成分是否可以把樣本比較不錯(cuò)的分開。

** 獲得訓(xùn)練數(shù)據(jù)前4個(gè)主成份的值 **

預(yù)測(cè)主成份的值，這里用的就是訓(xùn)練數(shù)據(jù)，所以得出訓(xùn)練數(shù)據(jù)主成分的值。

4 實(shí)戰(zhàn)二

R中自帶數(shù)據(jù)集data(Harman23.cor)數(shù)據(jù)集中包含305名受試者的8個(gè)身體測(cè)量指標(biāo)

5 進(jìn)階的主成分分析-psych包

正文中的princomp()函數(shù)為基礎(chǔ)包中進(jìn)行主成分分析的函數(shù)。 另外，R中psych包中提供了一些更加豐富有用的函數(shù)，這里列出幾個(gè)相關(guān)度較高的函數(shù)，以供讀者了解。

還有很多主成分分析結(jié)果可視化包，在直播我的基因組里面都提到過。

6 推薦一個(gè)R包factoextra

factoextra是一個(gè)R包，易于提取和可視化探索性多變量數(shù)據(jù)分析的輸出，包括：

• 主成分分析（PCA），用于通過在不丟失重要信息的情況下降低數(shù)據(jù)的維度來總結(jié)連續(xù)（即定量）多變量數(shù)據(jù)中包含的信息。

• 對(duì)應(yīng)分析（CA）是適用于分析由兩個(gè)定性變量（或分類數(shù)據(jù)）形成的大型應(yīng)變表的主成分分析的擴(kuò)展。

• 多重對(duì)應(yīng)分析（MCA），它是CA對(duì)包含兩個(gè)以上分類變量的數(shù)據(jù)表的適應(yīng)。

• 多因素分析（MFA）專門用于數(shù)據(jù)集，其中變量被組織成組（定性和/或定量變量）。

• 分層多因素分析（HMFA）：在將數(shù)據(jù)組織成層次結(jié)構(gòu)的情況下，MFA的擴(kuò)展。

• 混合數(shù)據(jù)因子分析（FAMD）是MFA的一個(gè)特例，專門用于分析包含定量和定性變量的數(shù)據(jù)集。

有許多R包實(shí)現(xiàn)主要組件方法。這些包包括：FactoMineR，ade4，stats，ca，MASS和ExPosition。然而，根據(jù)使用的包，結(jié)果呈現(xiàn)不同。為了幫助解釋和多變量分析的可視化（如聚類分析和維數(shù)降低分析），所以作者開發(fā)了一個(gè)名為factoextra的易于使用的R包。

7.主成分分析的生物信息學(xué)應(yīng)用

比如對(duì)千人基因組計(jì)劃的對(duì)VCF突變數(shù)據(jù)進(jìn)行主成分分析來區(qū)分人種：https://www./p/185869/

再比如每次做表達(dá)矩陣都需要根據(jù)樣本分組信息PCA分析及可視化看看分組是否可靠。

詳細(xì)例子見：http://github.com/jmzeng1314/GEO/ 

然后就是單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)也經(jīng)常會(huì)PCA看看分群，或者PCA來去除前幾個(gè)主成分因素來抹掉某些影響等等。

8. 主成分分析的其它可視化方法

看到一個(gè)包ropls可視化做的不錯(cuò)，本來以為ropls肯定是一個(gè)正常的r包，應(yīng)該是在cran里面，結(jié)果

現(xiàn)在什么包都往bioconductor里面丟真奇怪。但是作圖顏值還可以，大家可以看看：

后來仔細(xì)看標(biāo)題，終于明白了。

ropls: PCA, PLS(-DA) and OPLS(-DA) for multivariate analysis and feature selection of omics data

構(gòu)建的就是組學(xué)數(shù)據(jù)，所以放在bioconductor也是正常。

9.參考資料：

• http://www.cs./cosc453/student_tutorials/principal_components.pdf

• https://www.cs./picasso/mats/PCA-Tutorial-Intuition_jp.pdf

• http://www.cs./~samir/498/PCA.pdf

• http://www./ceo/Documentation/PCA_Outline.pdf

• http://people./~alawing/materials/ESSM689/pca.pdf （R相關(guān)）

• http://www2.dc./~cesar.souza/publications/pca-tutorial.pdf （2012）