Home | About | Posts | English

Principal component analysis (PCA) 的進階實作

之前的PCA實作以及主成分分析都建議先看一下，這篇文章則以直接上手操作可作為論文使用的功能為主。

資料輸入

將檔案存檔成csv檔之後放在你喜歡的路徑，依照下面的方式先設定工作路徑 (working directory) 再匯入資料，記得要將路徑的斜線設定正確。

setwd("D:/University/你喜歡的路徑")
dat <- read.csv("GeoConference_R.csv")

停！你的資料長怎樣

因為PCA的函數prcomp說笨不笨，說聰明也沒有到很聰明，你必須要清楚這個函數可以做什麼又不能作什麼，這篇文章是很好的教學材料，詳細解釋了許多原理，有想要深入研究的可以去讀。
但簡單來說，prcomp只能處理數字，而且絕對要數字，文字是無法處理的，如果你的資料長得像這樣，那麼你要進行PCA的可能就是後面的Al及Ca？總之了解你的目的即可。

Label	Horizon	Pedon	Al	Ca
SL1_A	A	1	30438	667
SL1_Bw	Bw	1	47903	525

進行PCA

依據你的資料有沒有文字會有不同，如果資料全部都是數字，那麼只要使用下面這段即可 (建議怕麻煩的人可以存成一個csv檔讀進來即可，簡單又快速)。

pca <- prcomp(dat, scale = TRUE)  
# scale = TRUE代表經過標準化，不用怕資料的變異不同！

但若你想把含有文字的檔案拿進來，之後還有要使用的話，由於prcomp中的參數formula需要為~ Al+Ca+Fe+P+K+Si+Ti+Cr+Pb+Mn這種形式，要將你想做PCA的數據title用加號連接起來，如果只有幾個的話用打字就好，如果有很多元素的話推薦使用下面這樣的程式碼。

paste(colnames(dat),collapse = "+")
# 會將你的資料title取出並以+連接
# 接著將上式的結果丟到下式取代這串Al+Ca+Fe+P+K+Si+Ti+Cr+Pb+Mn
pca <- prcomp(formula = ~ Al+Ca+Fe+P+K+Si+Ti+Cr+Pb+Mn,
              data = dat, 
              scale = TRUE)

那如果你的檔案很大、變數很多，就把檔案只留下要輸入的資料，然後丟下面這個程式碼。但要注意把方程式輸入成~.的形式就會自動將所有變數丟進去，這樣的缺點就是你只要有一個變數是文字或是包含缺值 (NA)，就會出現error喔!

pca <- prcomp(formula = ~.,
              data = dat, 
              scale = TRUE)

判讀資料與取得數據

接著直接進入我會看的東西，可以自己找需要的東西
Eigenvalue與累積比例

vars <- (pca$sdev)^2                # 提取特徵值
vars                                # 得到eigenvalue

props <- vars/sum(vars)             # 提取解釋百分比
cumulative.props <-  cumsum(props)  # 計算累積百分比
cumulative.props                    # 得到累積百分比

變數負荷量 (loading) 與樣本得分 (score)

write.csv(pca$rotation, "loading.csv") # 取得變數負荷
write.csv(pca$x, "score.csv")          # 取得樣本得分

存檔後的資料會在你原先設定的路徑裡面。

繪圖

由於我之前在PCA實作裡面寫的是使用ggfortify套件快速給予展示，可能不足以使用在論文中，底下分享我使用的繪圖程式碼，但可能比較複雜，拿前面的loading及score應該就足夠拿到SigmaPlot之類的軟體使用。我的示範資料大致長這樣：

Type	Al	Ca	Fe	P
1	265	786	30438	667
2	489	593	47903	525

在這個示範資料中，我有使用一個Type，這個資料並不是要用在PCA的資料，而是單純用於分類土壤，共分成四類，所以Type裡面有1、2、3、4，要讓R聽懂我們這四種分類的依據是將其定義為factor，使用下列程式碼：

dat$Type <- factor(dat$Type, levels = c("1", "2", "3", "4"))

這個level可以照自己設定的不一定要數字，例如化育層A、E、B、C也可以

autoplot(pca, data = dat, colour = 'Type', size=3.5, 
         loadings=F,            # 改成T可以顯示loading的箭頭
         loadings.label=F,      # 改成T可以顯示loading的名稱
         x=1, y=2,              # 修改x及y可以改XY軸是第幾個主成分
         scale=F)+              # 因為autoplot會自動改scale，而我們不想讓他改
  guides(colour=guide_legend(title="Type"))+
  scale_color_manual(labels=c("Type I", "Type II", "Type III", "Type IV"),   # 可以輸入你喜歡的名稱，化育層就輸入化育層
                     values=c("#F8766D", "#7CAE00", "#00BCF4", "#C77CFF"))   # 可以輸入你喜歡的顏色，這四個是我喜歡的顏色
                     # 但注意幾個分類這裡就要幾個，不然會錯誤

以上，祝大家都能得到自己喜歡的PCA結果！
備註：autoplot是從ggplot2延伸出來的函數，如果想要更加自由的操控可以直接使用ggplot2。