|
主成分分析（principal component analysis，PCA）是一種各個學科領域都廣泛使用的技術(shù)。也許很多讀者已經(jīng)學習過有關PCA的知識，網(wǎng)上也有很多講解PCA的文章，但這些文章的側(cè)重點不同，例如側(cè)重于線性代數(shù)、機器學習、生信應用等方面。然而，如何準確展示和解釋PCA的結(jié)果遠不是看上去那么簡單。Nature biotechnology 和 Nature reviews methods primers上都曾經(jīng)發(fā)文教大家相關的知識[1,2]。本期教程，我們將結(jié)合理論知識和作圖展示，教大家如何更好地掌握PCA。主成分分析（principal component analysis，PCA）是一種多變量統(tǒng)計方法，通過正交變換將一組可能存在相關性的變量轉(zhuǎn)換為一組線性不相關的變量，轉(zhuǎn)換后的這組變量叫主成分（principal components，PC）。 1.2 為什么要進行PCA呢？ 在許多領域的研究與應用中，往往需要對反映事物的多個變量進行大量的觀測。多變量大樣本無疑會為研究和應用提供了豐富的信息，但在多數(shù)情況下，許多變量之間可能存在相關性，不僅增加了問題分析的復雜性，也對分析建模帶來許多不便[3]。如果分別對每個指標進行分析，往往不能對問題產(chǎn)生綜合認識，而篩選并減少指標會損失很多信息，可能導致錯誤的結(jié)論。因此，需要找到一個合理的方法，在減少需要分析的指標同時盡量保留原指標中包含信息。由于各變量間存在一定的相關關系，那么我們就可以通過一種方式綜合原來的變量，建立少量的綜合指標來替代存在于各變量中的信息。上述過程是PCA實現(xiàn)數(shù)據(jù)降維的思路，從本質(zhì)上講，PCA旨在找到原始變量的線性組合（主成分），這些主成分可以解釋數(shù)據(jù)集中最大可能的變異量。 
圖1. 子刊上對于主成分分析的介紹。 文章“Principal component analysis”[1]中強調(diào)了發(fā)表文章需要進行的PCA工作流程，但這不是全部的細節(jié)。
圖2. PCA工作流程 其中大篇幅強調(diào)了雙標圖的制作與美化，這也是本教程想要強調(diào)的。我們先從PCA基本的要點講起。 （1）變量標準化 進行PCA的第一步就是變量標準化。這一步的目的是標準化連續(xù)初始變量的范圍，使得每個變量對分析的貢獻相等。因為如果初始變量的范圍之間存在很大差異，則范圍較大的變量將支配范圍較小的變量，這會導致結(jié)果有偏差，將數(shù)據(jù)轉(zhuǎn)換為可進行比較的尺度則可解決該問題。 通?？梢允褂脄-score標準化，即(數(shù)據(jù)-均值）/標準差（圖3）。z-score標準化可以將數(shù)據(jù)轉(zhuǎn)換為均值為0、方差為1的數(shù)據(jù)。 圖3. z分數(shù)標準化 此外，有時候一些R包會提到對數(shù)據(jù)進行中心化。中心化實際是數(shù)據(jù)-均值。這一步可將數(shù)據(jù)轉(zhuǎn)換為均值為0的數(shù)據(jù)。因此，可以把標準化看作是中心化的升級版。 在R語言中，這兩個操作都非常簡單，只需要使用R的基本函數(shù)scale()。 #構(gòu)造一個2列5行的矩陣，值為1到10x <- matrix(1:10, ncol = 2)#中心化scale(x, center = TRUE, scale = FALSE)#z-score標準化scale(x, center = TRUE, scale = TRUE)#或者直接使用默認參數(shù)執(zhí)行z-score標準化，和上一步結(jié)果一樣scale(x)
大多數(shù)用于PCA分析的R包都默認執(zhí)行中心化，例如FactoMineR包，只讓我們選擇是否標準化，在它們的函數(shù)里設置scale=TRUE就可以執(zhí)行標準化了，而無需我們提前算好。
降維過程是PCA的核心過程，需要計算特征向量（eigenvector）和特征值（eigenvalue）。特征向量與特征值總是成對出現(xiàn)，這是線性代數(shù)的概念，這里不再贅述。 “對于一個給定的線性變換A，它的特征向量經(jīng)過線性變換作用，得到的向量與原來的v保持同一直線， Av= λv，其中，λ是特征向量的長度在該線性變換下的縮放比例，是特征值。” 有兩種方法可以計算它們，一種是利用協(xié)方差矩陣進行特征值分解（eigenvalue decomposition，EVD），另一種是奇異值分解（singular value decomposition，SVD）?，F(xiàn)在流行的是SVD，因為它不要求數(shù)據(jù)是方陣，并且求解過程更加便捷。參考資料中對其進行了詳細介紹[1]，還畫了一個流程圖，這里就不展開講解了。圖3. PCA工作流程的框架圖。- 特征向量實際上是方差最大（信息最多）的軸的方向，稱之為主成分。
- 特征值與特征向量均為矩陣分解的結(jié)果。特征值表示標量部分，一般為某個主成分的方差，其相對比例可理解為方差解釋度或貢獻度 ；特征值從第一主成分會逐漸減小。
- 載荷（loading）是特征向量乘以特征值的平方根。主成分是所有原始變量的線性組合，載荷就是每個主成分上各個原始變量的權(quán)重系數(shù)。如果我們有三個原始變量X1、X2和X3，它們在第一個主成分中的載荷系數(shù)分別為0.5、0.3和0.8，則第一個主成分可以表示為0.5X1+0.3X2+0.8X3。
這幾個指標是我們進行圖表展示的基礎數(shù)據(jù)。計算過程就不展示了，因為我們可以使用R語言輕松算出它們。 在R語言中PCA常用的包及函數(shù)為： - stats 默認加載包，提供prcomp()與princomp()函數(shù)。
- psych 心理學研究包，提供principal()。
- FactoMineR 多元分析包，提供函數(shù)PCA()。
注意：函數(shù) princomp()使用特征值分解方法，而函數(shù) prcomp()和 PCA()使用奇異值分解的方法[4]，不同包的結(jié)果可能會有些微差別。但比較新的包都決定采用奇異值分解的方法。 圖4. prcomp和princomp函數(shù)的對比。 如果你閱讀過《R語言實戰(zhàn)》、《數(shù)量生態(tài)學》（《Numerical ecology with R》）之類的書籍，那么完成PCA分析并將結(jié)果可視化對你來說不是難事，你知道應該做哪些圖。這里我找到一個基于python的教程[5]，里面總結(jié)了常見的圖片。圖5. PCA分析的常用圖形。 在本教程中，我們會涉及上面的大多數(shù)圖片，但不是全部。我們的可視化方式也是基于R語言，其中一些圖片會有另外的展示方法。此外，當前環(huán)境與生態(tài)學領域非常流行一個R包叫做factoextraR，大多數(shù)PCA函數(shù)的結(jié)果都可以使用factoextraR包中提供的函數(shù)來輕松提取和可視化，我們后面會詳細介紹。
作為出版級文章的一部分，PCA分析最常用的圖可能是雙標圖（biplot）。這部分內(nèi)容涉及很多知識，下來我結(jié)合SAS官網(wǎng)上的博客文章[6]嘗試和大家再講一講。首先，當前兩個主成分解釋數(shù)據(jù)中大部分方差時，你可以通過將觀測值投影到前兩個主成分的范圍上來可視化數(shù)據(jù)。在PCA中，此圖稱為分數(shù)圖（得分圖）。你還可以將變量向量投影到主成分的范圍上，這稱為載荷圖。如下所示，PCA得分圖是一個散點圖，是用于可視化不同樣本之間相對位置的圖。在得分圖中，每個樣本被表示為一個點。通常，得分圖的橫坐標和縱坐標分別代表前兩個主成分的得分。PCA中，第一個主成分解釋了數(shù)據(jù)變異性的最大部分，第二個主成分解釋了第二大的變異性，以此類推。因此，得分圖的橫坐標和縱坐標所代表的數(shù)值大小是與數(shù)據(jù)集中的變異性和方差有關的。圖6. PCA分數(shù)圖。 載荷圖是幾個向量構(gòu)成的圖，顯示了變量和主成分之間的關系。如下所示，我們可以觀察到PetalWidth和PetalLength在x軸的投影較大，而在y軸上投影較小，說明它們對第一個主成分軸的貢獻高，而PetalWidth和PetalLength的夾角較小，說明它們密切相關。圖7. PCA載荷圖。 請注意，載荷圖的標尺（scale）比得分圖小得多。如果疊加這些圖，則矢量將顯示相對較小。在單個圖形中疊加得分圖和載荷圖的表達即為雙標圖，它由Gabriel (1971)首次在文獻中定義[7]。所有特征的觀測值將被投影到兩個維度上，從而使觀測值之間的距離近似保留（壓縮至兩個維度的主成分）。箭頭之間角度的余弦近似于變量之間的相關性，箭頭的長短近似于表示該變量的方差。下面顯示了一個示例。點是投影觀測值，向量是投影變量。
通常來說，它的左下軸是得分圖的刻度，上右軸是載荷圖的刻度。換句話說，左軸和底部軸是PCA得分圖的 — 使用它們來讀取樣本（點）的 PCA 分數(shù)。頂軸和右軸屬于載荷圖 — 使用它們來讀取每個特征（矢量）對主分量的影響程度。
下圖展示了如何閱讀一個雙標圖的信息[8]。
一個矩陣中，點是行（案例/樣本），向量是列（變量）。 兩個案例之間的距離近似于它們的相似性。 向量長度近似于變量的標準差。 兩個向量之間的余弦近似于變量之間的相關性。 案例在變量軸上的投影近似于最大值。
圖10. 雙標圖的解釋。 前面我們提到，得分圖和載荷圖通常具有不同的刻度。因此，在疊加得分圖和載荷圖時，需要重新縮放向量或觀測值（或兩者）。 根據(jù)分解過程中對列或行的保留，可以劃分行度量保留 （RMP） 雙標圖或列度量保留 （CMP） 雙標圖。在學術(shù)文獻中，這兩種類型分別稱為 JK 雙標圖和 GH 雙標圖。其中 JK 雙標圖更強調(diào)行的表示，而 GH 雙標圖更強調(diào)列的表示。為了生成對稱的雙標圖，我們需要平衡列和行的保留值，這就是所謂的SQRT雙標圖（也稱SYM雙標圖）。 無論你使用R、SAS、STATA還是別的什么軟件，雙標圖一般都會使用某種縮放，這個縮放方式和奇異值分解過程中的一個標量c有關，其取值在0~1之間，但不同軟件的計算方法可能不同，這個僅供大家了解，詳見參考資料[4-6]。 縮放有4種常見的選擇。每種縮放都側(cè)重強調(diào)觀測值對之間（如距離）、變量對之間（如角度）或觀測值與變量之間的某些幾何關系，而沒被強調(diào)的關系則會失真。由于所表示的數(shù)據(jù)通常具有高于二維或三維的維度（即其數(shù)據(jù)矩陣的秩），因此出現(xiàn)了如何在較低維度的空間中最佳地表示數(shù)據(jù)的問題。又稱“變量雙標圖”，是標量c=0時的圖。特點：向量被忠實地表示；觀測點之間的距離失真，不是歐氏距離；每個向量的長度與相應變量的方差成正比。這種雙標圖適合分析變量及其之間的關系。 在選擇保留變量關系時，觀測值可能投影到原點附近的微小區(qū)域，如下圖。
圖12. 變量雙標圖。 ② COV雙標圖 (Covariance biplot)GH雙標圖有一個改進版叫做協(xié)方差雙標圖，如圖8所示。它在標量c=0的基礎上還做了別的約束，使得每個向量的長度等于相應變量的標準差。這樣展示的效果比GH雙標圖更好[6]。
特點：變量的標準差和相關性被較好表達。如果變量已預先標準化為標準差等于1，則通常會在COV雙標圖上繪制一個單位圓，因為變量點的長度都小于或等于1——變量點離單位圓越近，顯示得越好。 又稱“距離雙標圖”或者“形式雙標圖（form biplot）”，是標量c=1時的圖。特點：觀測值被忠實地表示，觀測值之間的距離是歐氏距離；向量之間的角度會因縮放而失真。這種圖和PCA得分圖一樣關注樣本之間的關系，是最為常用的類型。也就是，這種雙標圖適合分析觀測值及其之間的關系。
圖13. 距離雙標圖。 JK雙標圖和COV雙標圖被稱為對偶雙標圖（dual biplots）：每一個都是另一個的對偶[9]。這是很多軟件中常用的2種選擇。由于每個主成分的系數(shù)平方和為1。有時候，為了在圖形中表達這一點，軟件會將變量的線長“拉伸”到其原始長度（即使用變量對應的載荷系數(shù)作為其坐標）。圖14. SYM雙標圖。 ④ SYM雙標圖（symmetric biplot）又稱“對稱雙標圖”，是標量c=1/2時的圖。特點：折衷處理觀測值和向量，兩者都一定程度失真。這種雙標圖適合分析觀測點和變量之間的關系。似乎很少有人使用這種方法。
圖15. SYM雙標圖。 以下是我從網(wǎng)上找的不完全的總結(jié)。 圖16. 3種類型的雙標圖。 此外，《數(shù)量生態(tài)學》一書中使用vegan包的rda()函數(shù)進行主成分分析，書中提到常用的2種類型的縮放，或者說標尺，分別對應距離雙標圖和相關雙標圖。下面的框中展示了具體細節(jié)，但這些只是在vegan包中的方法。 按照vegan包的邏輯，如果主要目的是解釋對象之間的關系，則可將PCA排序圖執(zhí)行I型縮放，即I型標尺（scaling 1）-距離雙標圖（distance bioplot）。特征向量被標準化為單位長度，關注的是對象之間的關系。雙標圖中對象之間的距離近似于多維空間內(nèi)的歐式距離，代表變量箭頭之間的角度沒有意義。 變量之間的相關性是通過變量向量之間的角度，而非向量頂點之間的距離來描述的。若期望關注變量間的相關性，需將PCA排序圖執(zhí)行II型縮放，即II型標尺（scaling 2）-相關雙標圖（correlation bioplot）。特征向量被標準化為特征根的平方根，關注的是變量之間的關系。雙標圖中對象之間的距離不再近似于多維空間內(nèi)的歐式距離，代表變量箭頭之間的夾角反映變量之間的相關性。例如下面這個圖。案例分析了將探討濕地植物不同群落的環(huán)境變量（主要是水化學）的相互關系。
圖17. 兩種標尺的雙標圖。左圖使用標尺1（關注樣本之間的距離），右圖使用標尺2（關注物種/變量之間的相關性，反映在特定向量的角度）。左圖中的圓是所謂的平衡貢獻圓。 順便一提，在距離雙標圖中，可以通過平衡貢獻圈（circle of equilibrium contribution）評估變量的相對重要性（如上面左圖所示）。
平衡貢獻圈半徑為(d/p)的平方根，其中d是所展示的PCA軸數(shù)（通常d = 2），p是變量個數(shù)，可知半徑長度是一開始便定義好的固有特征。平衡貢獻圈的半徑代表變量的向量長度對排序的平均貢獻率。如果某個變量的箭頭長度大于圓的半徑，代表它對這個排序空間的貢獻大于所有變量的平均貢獻，可以認為這些變量是相對更重要的變量。 vegan包的官方說明中提到，包中的縮放或RDA結(jié)果與大多數(shù)其他軟件包不同，由于其細節(jié)非常復雜，讀者可以從文檔“Design decision and implementation details in vegan”[10]查看相關的解釋。提醒大家，在不清楚一個R包的縮放邏輯時，一定要先看幫助文檔，以免誤讀PCA的結(jié)果。以上是一些重要的知識，其他知識我們將在實戰(zhàn)中介紹。還是來看R語言的實戰(zhàn)例子吧~ 2.1 前期準備 這里我們使用FactoMineR包來做PCA，并使用factoextra來提取并可視化結(jié)果。首先我們要安裝這兩個包。
#安裝R包install.packages("factoextra")install.packages("FactoMineR")
#加載包library(FactoMineR)library(factoextra)
這里我們用的是數(shù)據(jù)是R自帶的iris數(shù)據(jù)集，在以前的教程中我經(jīng)常用它做示范。該數(shù)據(jù)包含3個鳶尾花物種的50個樣本，且每個樣本中測量了4個特征，即萼片和花瓣的長度和寬度。
圖18. iris數(shù)據(jù)集的信息 現(xiàn)在我們打算用這4個特征來做PCA。 #加載數(shù)據(jù)data(iris)#查看數(shù)據(jù)結(jié)構(gòu)str(iris)
PCA的第一步是數(shù)據(jù)標準化。在FactoMineR的PCA()函數(shù)中，默認設置對數(shù)據(jù)執(zhí)行標準化，即scale.unit=TRUE。該函數(shù)用法如下：
我們可以直接進行PCA。結(jié)果會返回我們2張不同的圖。
#PCA只對數(shù)值矩陣生效，我們不需要第5列（物種）res.pca <- PCA(iris[,-5])
上面2張圖分別為PCA的得分圖和載荷圖（變量圖）。可以看到，第一個PC解釋了72.96%的方差，第二個PC解釋了22.85%的方差。 而具體的結(jié)果我們可以通過summary()得到。
可以使用factoextra包的get_pca_var()函數(shù)來提取變量的分析結(jié)果。 # 提取變量的分析結(jié)果var <- get_pca_var(res.pca)var
coord表示用于創(chuàng)建散點圖的變量坐標。coord實際上就是成分載荷，指觀測變量與主成分的相關系數(shù) cor表示相關系數(shù) cos2表示因子質(zhì)量，其計算公式為：var.cos2 = var.coord * var.coord contrib表示包含變量對主成分的貢獻（百分比），其計算公式為： (var.cos2 * 100) / (total cos2 of the component)。
此外，我們可以使用factoextra包別的函數(shù)來有針對性地提取別的結(jié)果，基本函數(shù)有這些： get_eigenvalue(res.pca)：提取主成分的特征值/方差 fviz_eig(res.pca)：可視化特征值 get_pca_ind(res.pca)，get_pca_var(res.pca)：分別提取個體和變量的結(jié)果。 fviz_pca_ind(res.pca)，fviz_pca_var(res.pca)：分別可視化結(jié)果個體和變量。 fviz_pca_biplot(res.pca)：制作主成分分析散點圖biplot圖。
為了更直觀地展示每個主成分解釋的方差，我們可以使用函數(shù)fviz_screeplot()繪制所謂的碎石圖(screeplot)。這個函數(shù)有一些美化參數(shù)，大家自己去看看。
#繪制主成分碎石圖，查看每一個主成分能在多大程度上代表原來的特征fviz_screeplot(res.pca, addlabels = TRUE)
隨之而來的問題是，主成分的個數(shù)如何確定？簡單而言，特征值和累積貢獻率，都是判斷提取幾個主成分的標準。一方面，我們通常提取特征值大于1的因素（這稱為Kaiser-Guttman 準則，如果一個因子的特征值小于1，說明其比原始變量解釋的方差還少）；另一方面，我們希望提取的主成分累積解釋超過80%的方差。假如前兩個PC能解釋超過80%的方差（經(jīng)驗標準），這時候我們就能很有底氣地說提取前2個PC就行了。但這不是死標準，生態(tài)學中沒有被廣泛接受的客觀方法來決定有多少主成分就足夠了。一行簡單的代碼可以幫我們提取特征值。只需要使用get_eigenvalue()函數(shù)。
我們也可以把特征值做成碎石圖。 #繪制主成分碎石圖(以特征值為縱軸)fviz_screeplot(res.pca, choice = "eigenvalue",addlabels = TRUE)
按這么看來，我們應該提取前2個主成分就行了（實際上我們自己的數(shù)據(jù)沒這么好的效果）。
我們還可以通過碎石圖（以解釋方差比例為縱軸）來判斷“拐點”，即觀察從哪個主成分后下降的趨勢就趨于平穩(wěn)了，換句話說，當主成分數(shù)目從k增加到k+1時，如果特征值出現(xiàn)較劇烈的下降，即第k+1個點為拐點，表示相比較于前k個主成分所提取的變量信息，第k+1個主成分只包含了非常少量的變量信息，因此提取k個主成分較為合適。當然，這種觀察也比較主觀。于是，有人提出了斷棍模型（broken stick model）來輔助我們判斷。斷棍模型的主要原理是將單位長度的棍子隨機分成與PCA軸數(shù)一樣多的幾段，然后將這些斷棍按照長短依次賦予對應的軸（即最長的棍子賦予第一軸，第二長的賦予第二軸，依此類推）。公式如下[11]：
選擇標準的標準有兩種，一是選取特征根大于所對應的斷棍長度的軸，第二種是選取特征根的總和大于所對應斷棍長度總和前幾軸。R語言自帶的stats包的函數(shù)princomp()和prcomp()生成的PCA結(jié)果可以使用函數(shù)screeplot()來展示斷棍模型的圖。而factoextra包沒有開發(fā)這個功能（可能是覺得不太重要）。我們這里需要自己編寫下繪制斷棍模型的代碼。
#提取特征根res <- res.pca$eig[,1]#斷棍模型n <- length(res)bsm <- data.frame(j=seq(1:n), p=0)bsm$p[1] <- 1/nfor (i in 2:n) { bsm$p[i] = bsm$p[i-1] + (1/(n + 1 - i))}bsm$p <- 100*bsm$p/nbsm#可視化barplot(t(cbind(100*res/sum(res),bsm$p[n:1])), beside=TRUE, main="% 變差", col=c("bisque",2), las=2)legend("topright", c("% 特征根", "斷棍模型"), pch=15, col=c("bisque",2), bty="n")
這么看來，選擇1個主成分都行，選2個更行。補充一句，有時候，第二個主成分和第三個主成分差距不大，我們也可以用PC1和PC3來作圖。正如前面說的，選擇幾個主成分沒有一個統(tǒng)一的標準，需要結(jié)合多個準則并根據(jù)你具體探究的問題來決定。2.4 PCA結(jié)果的得分圖 前面我們看到了，函數(shù)PCA()會返回我們2張圖，即PCA得分圖/觀測值坐標圖和載荷圖。但這是作為探索性分析默認的顯示方式。當我們確定了要提取前兩個主成分，我們就得好好畫下這個圖。 觀測值坐標圖對應了factoextra包中的fviz_pca_ind()函數(shù)，而變量分布圖對應fviz_pca_var()函數(shù)。 首先，我們畫得分圖。這里我們使用物種信息來著色和分組。 fviz_pca_ind(res.pca, geom.ind = "point", # 只顯示點而不顯示文本，默認都顯示 # 設定分類種類 col.ind = iris$Species, # 設定顏色 palette = c("#00AFBB", "#E7B800", "#FC4E07"), addEllipses = TRUE,# 添加置信橢圓 legend.title = "Groups",)
置信橢圓定義的區(qū)域包含了95%的樣本。假如我們想看全部樣本的邊界，可以把它改成多邊形。 # PCA得分圖fviz_pca_ind(res.pca, mean.point = F, # 去除分組的中心點，否則每個群中間會有一個比較大的點 label = "none", # 隱藏每一個樣本的標簽 habillage = iris$Species, # 根據(jù)樣本類型來著色 palette = c("purple","orange","blue"), # 三個組設置三種顏色 addEllipses = TRUE, # 添加邊界線 ellipse.type = "convex" # 設置邊界線為多邊形)
有的人可能想看每一個樣本的標簽，那么我們可以不隱藏文本。
fviz_pca_ind(res.pca, # 設定分類種類 col.ind = iris$Species, # 設定顏色 palette = c("#00AFBB", "#E7B800", "#FC4E07"), addEllipses = TRUE,# 添加置信橢圓 legend.title = "Groups",)
上圖得調(diào)下文字大小和位置，不然太亂了。 還可以根據(jù)對樣本的質(zhì)量/貢獻改變其大小。 fviz_pca_ind(res.pca, pointsize = "cos2", pointshape = 21, fill = "#E7B800", repel = TRUE # 避免文本相互重疊)
另外就是繪制變量分布圖。最簡單的圖只需要一行代碼。
在factoextra包中，這個圖被稱為相關圓（Correlation circle）。factoextra包的作者kassambara在官方教程里寫道[12]，變量和主成分（PC）之間的相關性被用作PC上變量的坐標。注意：變量的表示與觀測值的圖不同。觀測值由其投影表示，但變量由其相關性。 因此，在相關圓上，不會有變量的長度超過1。 相關圓圖顯示了所有變量之間的關系，其要點是： 正相關變量之間夾角小于90°，而負相關變量之間的夾角大于90°，夾角等于90°的變量之間沒有相關性（實際上接近90°的就可以認為不太相關）。 負相關變量位于圖原點相對的兩側(cè)（相對的象限）。 變量與原點之間的距離用于測量圖上變量的質(zhì)量（表示為cos2）。遠離原點的變量在圖上得到了很好的表示。
fviz_pca_var(res.pca, col.var = "contrib", #根據(jù)貢獻度著色 gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07") #設置色帶)
我們可以看主成分分析結(jié)果中的contrib的數(shù)值高低來觀測變量對主成分的貢獻，更直觀的方法是做熱圖。
# 獲取樣本的主成分分析結(jié)果var <- get_pca_var(res.pca)#加載繪制熱圖的包library("corrplot")#繪制熱圖corrplot(var$contrib, is.corr = FALSE)
類似的，我們可以使用cos2這個指標，cos2反映了各個主成分中各個變量的代表性。一個變量的所有主成分cos2值加起來等于1。對于主成分而言，某個變量的cos2越接近1，則說明變量對該主成分的代表性越高；cos2越接近0，則說明變量對該主成分的代表性越差。 fviz_pca_var(res.pca, col.var = "cos2", #根據(jù)代表性著色 gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07") #設置色帶)
# 獲取樣本的主成分分析結(jié)果var <- get_pca_var(res.pca)#繪制熱圖corrplot(var$cos2, is.corr = FALSE)
我們也可以使用函數(shù)創(chuàng)建條形圖，來看變量cos2反映的代表性。 # Total cos2 of variables on Dim.1 and Dim.2fviz_cos2(res.pca, choice = "var", axes = 1:2)
較高的 cos2 表示變量在主成分上有良好的代表性。在這種情況下，變量的位置靠近相關圓的圓周。反之，低 cos2 表示變量未由 PC 完美表示。在這種情況下，變量靠近圓心。 那如果 2.6 PCA結(jié)果的雙標圖 要繪制一個簡單的雙標圖，只需要使用如下代碼： fviz_pca_biplot(res.pca, repel = TRUE, col.var = "#2E9FDF", # 變量顏色 col.ind = "#696969" # 樣本顏色)
先給大家講怎么美化圖片，后面我們再說這個圖使用的縮放類型。
這個圖的樣本標簽還是隱藏比較好。 # 只保留變量的標簽fviz_pca_biplot(res.pca, label ="var")
我們可以在此基礎上增加分組和置信橢圓。
fviz_pca_biplot(res.pca, col.ind = iris$Species, addEllipses = TRUE, label = "var", col.var = "black", repel = TRUE, legend.title = "Species")
從該圖我們可以看出，Petal.width和Petal.Length對于確定PC1評分和區(qū)分物種組很重要。由于該圖片是ggplot2對象，我們可以使用和ggplot2一樣的美化思路去美化它。比如，更改一下配色。fviz_pca_biplot(res.pca, col.ind = iris$Species, addEllipses = TRUE, label = "var", col.var = "black", repel = TRUE, legend.title = "Species")+ ggpubr::fill_palette("jco")+ #樣本填充顏色，使用jco期刊的風格 ggpubr::color_palette("npg") #變量顏色，使用npg期刊的風格
關于更多圖片美化的信息可以見參考資料[12]，這是factoextra包的作者kassambara寫的官方教程。現(xiàn)在的問題是，它使用了什么縮放？ 作者kassambara提醒我們： 僅當數(shù)據(jù)集中的變量和個體數(shù)量較少時，雙標圖才可能有用；否則最終的plot將是不可讀的。 另請注意，個體和變量的坐標不是在同一空間上構(gòu)造的。因此，在雙標圖中，您應該主要關注變量的方向，而不是它們在圖中的絕對位置。
粗略地說，雙標圖可以解釋如下： 其幫助文檔中并沒有明確說明其所用方法，也沒有設置靈活的參數(shù)來改變標尺。答案并不像我們想象中那么明了。但我們將該圖和之前的JK雙標圖對比，可以發(fā)現(xiàn)坐標刻度和樣本的位置是一致的，只是箭頭長度不同。既然作者說我們應該關注變量的方向，而不是絕對位置，那么這個圖應該是一種JK雙標圖的變體，其箭頭長度被等比放大了。為了美觀，一些別的軟件也會這樣做。 如果你只是想簡單地描述biplot圖，記住JK雙標圖更適合分析觀測值之間的關系就行。其他關系是有一定失真的，但是也能給我們提供大致的信息。此外，圖片的縱橫比會對雙標圖的角度有一定影響，大家可以拉伸下圖片看看變化。因此，如果我們有多組雙標圖要繪制，甚至比較，我建議大家固定一下輸出坐標系的縱橫比為1:1，使用原本的繪圖+ coord_fixed()。如果有讀者了解更多factoextra包中相關的知識，還請不吝賜教~這個包我們就介紹到這兒。用不同包做出來的雙標圖是否會差別很大呢？ 可以想象，箭頭長度的差異也許很大，因為JK雙標圖經(jīng)常被設置箭頭的拉伸。但有時你會發(fā)現(xiàn)不同包作圖結(jié)果中，變量角度上也有差異。例如，有人在stackexchange上提問“Why do arrows of PCA graph have different angles between biplot and ggplot functions?”[13]。提問者認為factoextra包的fviz_pca_biplot()和基本函數(shù)biplot()給出的圖角度差異很大。 這也許不是軟件出錯了。評論區(qū)的網(wǎng)友認為首先要排除不同包中數(shù)據(jù)標準化的方式一致，并且圖片的縱橫比是一致的。在排除這種差異后，提問者發(fā)現(xiàn)箭頭的角度沒啥差異了，但箭頭在樣本空間的位置有差異。
提問者提供的圖 這種差異的原因很可能是上面左右兩個圖都是兩個獨立圖的疊加（使用了不同的縮放方法）。我們可以看到，左圖的刻度和PCA得分圖并不一致，應該是縮放了。因此將箭頭與點進行比較時的任何差異都是由于對每個圖進行歸一化的不同約定造成的。即使是JK雙標圖，也會有人對得分縮放嗎？是的，因為PC1和PC2解釋的方差可能差異很大，數(shù)據(jù)范圍跨度很大，為了使其合適的展示在一個圖中，有的軟件（包）會設置某種縮放。我注意到stackoverflow上相關的提問“How to create a biplot with FactoMineR?”，評論區(qū)提到了標尺的問題。看來這是學習過biplot的人都知道的注意事項，不是那種無關緊要的知識。網(wǎng)上的一條評論 令人費解的是，標尺的參數(shù)在R的原始函數(shù)biplot.princomp()中是有定義的，幫助文檔里也有一些解釋，為什么很多后來開發(fā)的R包并沒有允許我們選擇GH、JK等標尺呢？
當然，如果實在沒有好用的包，我們也可以使用princomp()做PCA，用biplot()做雙標圖。但原始的函數(shù)作圖不是很方便，不信大家可以去看看我找到的某篇教程[14]。biplot函數(shù)作的圖是真的不好看吧？ 總之，可能這個問題之前真沒啥人關注（也許現(xiàn)在的人很少從PCA的圖片去解釋這些信息？）但準確選擇雙標圖的類型能幫助我們更好地推斷結(jié)論。3.3 探索更好的雙標圖 在潛心研究了網(wǎng)上的ggbiplot、ggfortify、ggord、PCAtools等若干個包之后，我發(fā)現(xiàn)這些包在某些方面很優(yōu)秀，并且它們都可以繪制雙標圖，但不一定好用，其標尺的選擇和設置也不一定很明了（大家可以幫我研究下是不是這樣）。 接下來，我將向大家展示有關雙標圖的更多內(nèi)容。 （1）ggord包作雙標圖 ggord包是一個專門用來繪制排序圖的包[15]。它基于ggplot2開發(fā)。需要執(zhí)行以下代碼來安裝它。 # Enable the r-universe repooptions(repos = c( fawda123 = 'https://fawda123.v', CRAN = 'https://cloud.'))
# Install ggordinstall.packages('ggord')
我閱讀了該包的誕生過程，一篇2015年的博客[15]，作者說他是“重新造輪子”，盡管已經(jīng)有ggbiplot和factoextra珠玉在前，這倆包已經(jīng)提供了幾乎完全覆蓋 R 中多變量和排序分析的繪圖結(jié)果。但作為一個固執(zhí)的人，他不愿放棄自己的包，并開始探索改進這些現(xiàn)有包中雙標圖的一些功能。處于主成分分析，他還將雙標圖拓展到用于nonmetric multidimensional scaling, multiple correspondence analysis, correspondence analysis 和 linear discriminant analysis等。 用prcomp()、PCA()等函數(shù)得到的PCA結(jié)果可以直接傳遞給它。也就是上一節(jié)中我們命名為res.pca的變量。我們可以使用以下代碼畫雙標圖。 p<- ggord(res.pca, iris$Species)p
為了顯示完整，我們可以調(diào)整下xlim和ylim參數(shù)來確定繪圖區(qū)間。p <-ggord(res.pca, iris$Species,xlim=c(-3,3.5),ylim=c(-3,3))p
這個圖看起來應該是JK雙標圖。它的箭頭長度較短。 函數(shù)ggord()內(nèi)置的參數(shù)非常多，就不逐一介紹了，大家可以看看幫助文檔。 似乎ggord不允許我們縮放向量的長度，但我們可以改變箭頭的一些屬性，包括其它的屬性，如點的大小，橢圓的透明度等。
p <- ggord(res.pca, iris$Species, arrow = 0.4,#箭頭頭部的長度 vec_ext = 1,#箭頭標注的遠近 veccol = 'red', #箭頭的顏色 veclsz = 1,#箭頭的粗細 vectyp = 'longdash',#箭頭的線型，這里改成長虛線 size = 3,#點大小 alpha_el=0.3)#置信橢圓透明度p
（2）ggbiplot包作雙標圖 ggbiplot包是一個比較老的包，很多功能沒更新了。library(devtools)install_github("vqv/ggbiplot")
使用同名函數(shù)ggbiplot()來繪制結(jié)果。
p<-ggbiplot(res.pca, groups = iris$Species, #設置分組變量，點的顏色會根據(jù)分組改變 ellipse = TRUE, #繪制置信橢圓 circle = TRUE) #繪制相關圓p
看得出來，得分軸被進行了標準化。
其中scale參數(shù)就是調(diào)整雙標圖的縮放方式。當scale=1，展示方式為協(xié)方差雙標圖，變量之間的內(nèi)積近似于協(xié)方差，點之間的距離近似于馬氏距離；當scale=0，展示方式為形式雙標圖。但是如果我們只改這個參數(shù)，另scale=0，會得到報錯。網(wǎng)上有個老外說作者可能放棄開發(fā)這個函數(shù)了。這人叫Kevin Blighe，他是一個新的R包PCAtools包的作者。（3）PCAtools包作雙標圖 PCAtools包提供通過PCA進行數(shù)據(jù)探索的功能，并允許用戶生成可供出版的圖形。用戶還可以使用它確定主成分的最佳數(shù)量、繪制不同主成分的成對圖等。此外，這個包似乎是個組學神器。 需要執(zhí)行以下代碼來從github中安裝它。 #安裝包if (!requireNamespace('devtools', quietly = TRUE)) install.packages('devtools')devtools::install_github('kevinblighe/PCAtools')#加載包library(PCAtools)
PCAtools包中的函數(shù)biplot()可以用于繪制雙標圖。首先，我們需要使用該包的函數(shù)pca()來得到主成分分析的結(jié)果，這是一個"pca"類型的對象，函數(shù)biplot()只能接受這樣的輸入對象。所以，這會兒大家就暫時先忘了前面那些包吧~
參數(shù)mat是用于PCA的輸入數(shù)據(jù)，而metedata是輔助數(shù)據(jù)。這兩個數(shù)據(jù)需要具有行列名，使得rownames（metadata）==colnames（mat）。值得注意的是，一定要設置transposed=TRUE，代表mat矩陣是轉(zhuǎn)置的。因為這個包默認mat的行是變量，列是樣本，這和我們平時保存數(shù)據(jù)的方式相反。我認為這里對于行列名的嚴格要求制造了很多麻煩。至今我在試用的時候還會出現(xiàn)很大出人意料的錯誤。但是這個包的其他功能看上去很酷！比如，繪制每一對主成分組成的矩陣圖，這個很適合探索性分析。核心函數(shù)是pca()。其中，scale參數(shù)需要設置成TRUE，意思是需要對數(shù)據(jù)使用標準化。這個FactoMineR包的PCA()函數(shù)和PCAtools包的pca()函數(shù)默認對數(shù)據(jù)進行中心化，所以這個設置我們可以不管（這里寫出來以便閱讀）。 #為數(shù)據(jù)集設置行名row.names(iris)<-1:150#執(zhí)行PCAres <- pca(iris[,-5], scale=TRUE, center=TRUE, transposed=TRUE)#查看變量的類型class(res)
輸出數(shù)據(jù)的方差。 #查看方差res$variance#也可以使用該包的函數(shù)getVars(res)
此外，還有一些簡單的函數(shù)幫助我們提取PCA的部分結(jié)果。 #提取主成分getComponents(res)#提取載荷getLoadings(res)
我們可以使用函數(shù)screeplot()來繪制碎石圖。
screeplot(res, axisLabSize = 18, titleLabSize = 22)
條形圖顯示了每個主成分解釋的方差占比。折線圖是累積解釋的方差占比。相應的數(shù)值是和我們前面用別的包做的PCA結(jié)果一樣的。
作者提到，人們常說的“biplot”存在不同的解釋。在OMIC時代，對于大多數(shù)普通用戶來說，雙標圖是二維空間中樣本的簡單表示，通常只關注前兩個PC。然而，Gabriel KR（Gabriel 1971）對雙標圖的最初定義是在同一空間中繪制變量和觀測（樣本）的圖。變量由從原點繪制的箭頭指示，箭頭指示它們在不同方向上的'weight’。使用biplot()繪制一個普通的樣本散點圖的代碼如下： 當我設置更多參數(shù)時，錯誤不斷。所以我暫時放棄了。（4）ggfottify包作雙標圖 最后，我找到了一個看起來維護很好的包（更新于2023-03-20），ggfortify包。據(jù)cran官網(wǎng)介紹，ggfottify基于ggplot2開發(fā)，是一個強大的繪圖百寶箱，可以作出常用統(tǒng)計繪圖，如GLM、 時間序列、PCA分析、聚類和生存分析的圖。但就雙標圖而言，它局限性還是挺明顯。只支持prcomp()和 princomp()函數(shù)的結(jié)果作為輸入。那么，我們使用prcomp()函數(shù)來得到PCA的結(jié)果。當然，這個結(jié)果和前面基本一致。res1 <- prcomp(x = iris[ , -5], center = TRUE, scale. = TRUE)
接下來，使用它的繪圖函數(shù)autoplot()，設置loadings=TRUE意思是顯示變量。顧名思義，這個函數(shù)的內(nèi)置參數(shù)應該比較少，它傾向于自動出圖，至少我在說明文檔里沒找到它其他參數(shù)的介紹。
#繪圖autoplot(res1, data = iris, colour = 'Species', loadings = TRUE)
這個圖和我們之前做的上下顛倒，是正常的，因為prcomp()函數(shù)得到的載荷系數(shù)和我們之前用的包剛好正負相反，這個不會影響別的什么。 但是這個是哪種雙標圖呢？考慮到其對prcomp()函數(shù)的依賴，這個可能是和biplot()的結(jié)果一樣的。 總的來說，我認為這幾個包中g(shù)gord在美化上更容易，可以畫出更好看的雙標圖。而PCAtools提供了更多的PCA分析工具，在數(shù)據(jù)探索方面很有前景，但雙標圖的設置太復雜而存在很多使用門檻。其他2個包不太適合用于PCA分析的出圖。 [1] Principal component analysis (https:///10.1038/s43586-022-00184-w)[2] What is principal component analysis? (https:///10.1038/nbt0308-303)[3] 主成分分析(PCA)原理分析&Python實現(xiàn) (https://blog.csdn.net/dongke1991/article/details/126774638)[4] Principal Component Analysis in R: prcomp vs princomp - Articles - STHDA (http://www./english/articles/31-principal-component-methods-in-r-practical-guide/118-principal-component-analysis-in-r-prcomp-vs-princomp/)[5] Example gallery — psynlig 0.2.1.dev0+1b3d658 documentation (https://psynlig./en/latest/auto_examples/index.html#principal-component-analysis)[6] What are biplots? - The DO Loop (https://blogs./content/iml/2019/11/06/what-are-biplots.html)[7] Gabriel, K. R. (1971). The biplot graphical display of matrices with applications to principal component analysis. Biometrika, 58, 453–467. doi:10.2307/2334381.[8] On the use of biplot analysis for multivariate bibliometric and scientific indicators (https:///10.1002/asi.22837)[9] Multivariate Statistics (https://www./microsite/multivariate-statistics/biplots.html)[10] Design decisions and implementation (https://cran./web/packages/vegan/vignettes/decision-vegan.pdf)[11] Principal Component Selection: The Broken-Stick Model | Neuro Mo (https://www./post/broken_stick/)[12] PCA - Principal Component Analysis Essentials - Articles - STHDA (http://www./english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/#eigenvalues-variances)[13] r - Why do arrows of PCA graph have different angles between biplot and ggplot functions? - Cross Validated (https://stats./questions/580490/why-do-arrows-of-pca-graph-have-different-angles-between-biplot-and-ggplot-funct)[14] Biplot using base graphic functions in R (http://agroninfotech./2020/06/biplot-for-pcs-using-base-graphic.html)[15] Reinventing the wheel for ordination biplots with ggplot2 | R-bloggers (https://www./2015/05/reinventing-the-wheel-for-ordination-biplots-with-ggplot2/)成體系地講解PCA確實不容易，不同學科的人對它可能有不同的期待。我只能盡量把重點講到。如果大家看完這個教程還有很多疑問，歡迎來交流群和我討論~其實在文獻中還有一些PCA的高級可視化方案，我們之后的教程中會講到。作為頂刊可視化的系列教程，本公眾號的后續(xù)推文將講介紹如何復現(xiàn)一篇Nature communications上的PCA圖。
如果覺得文章對您有用請點擊右下角“在看”，讓文章能夠幫助到更多的小伙伴。聲明：參考材料可能包含大量信息，如您對文章內(nèi)容很感興趣，建議您閱讀推文中列出的參考材料。本文涉及的部分譯文內(nèi)容僅供參考，一切以英文原文為準。由于編譯者水平有限，不當之處敬請諒解。
|
