作圖的目的是希望在圖里面發(fā)現問題或者解釋問題，當然更本質一點就是你想解決什么問題？

前幾天做了一個PCA的圖，圖是畫出來了，但是問題有很多，比如說主成分是是啥意思，圖里面的箭頭有什么含義？為了不做無意義的重復，所以寫一篇文章嘗試做一個解釋。

我們以R語言自帶的數據集iris作為例子來演示。

data(iris)

iris翻譯成中文就是鳶(yuan， 第一聲)尾花（如下圖）, 我建議你在R語言里用?iris了解更多這個數據集的出處。

iris

先讓我們用head簡單看下這個數據集的前面幾行。其中前面四列是一朵鳶尾花的一些形態(tài)特征衡量值，自行翻譯各個單詞的含義。最后一列是這朵花的所屬物種，根據分類，這些花來自于"setosa, versicolor,virginica"

> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

在平常人眼里, 看到這些數值可能沒有任何想法，也不知道這到底是什么品種，但是對于一個研究鳶尾花的人，這些數值可能會在他們的腦海里重構出一朵鳶尾花，然后迅速判斷出這是什么品種。我和閱讀文章的人都一樣，也不知道如何區(qū)分鳶尾花的品種，如果想去研究這種規(guī)律，應該怎么做呢？

人類的學習過程過程大多數多看多提出一些模型規(guī)律，然后基于這個規(guī)律去判斷，如果錯了改正模型，直到找到一個好用的標準為止。因此，后面我們也就是在當前的數據集下多試試，看看能不能提出一些模型。

可視化是非常強大的工具，讓我們先看看能不能用“ Sepal.Length”和“ Sepal.Width”進行區(qū)分

library(ggplot2)
ggplot(iris,aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(colour=Species))

fig1

從上圖中，我們發(fā)現根據萼片(sepal)的寬度和長度似乎能夠區(qū)分出"setosa"和"versicolor","virginica"，但是"versicolor","virginica"不太好分。

讓我們再試試花瓣(petal)的長度和寬度。從下圖，我們可以發(fā)現在這兩個屬性下，不同品種的鳶尾花似乎有了明顯的分界線。

ggplot(iris,aes(x=Petal.Length, y=Petal.Width)) + geom_point(aes(colour=Species))

image

我們可以建立一個模型就是

• 當花瓣長度小于2時就一定是setosa
• 當花瓣的寬度在0.75~1.75之間，花瓣的長度在3~5，則是versicolor,
• 當花瓣的寬度大于1.75時，花瓣的長度大于5，則是virginica

這個模型依舊還不完善，存在一些點難以區(qū)分。不過我們還可以以花瓣長度和萼片長度，花瓣寬度和萼片寬度等你能想到的組合進行作圖，說不定能找到更好的模型。更好的情況是，如果人類能夠想象出一個四維的空間，在這個空間里面同時展現這四個變量，或許就能更加直觀的找到規(guī)律，可惜，我們目前還沒有這種能力

有沒有一種方法能夠把這種多維數據降維，也就是抓住數據集的主要矛盾呢？當然是有，這個方法就是100多年前，由卡爾·皮爾森提出的 主成分分析，通過將原來眾多具有一定相關性的變量，重新組合成一組新的互相無關的綜合變量。比如說我們可以將花瓣長度和寬度組合成花瓣的面積，當然實際計算可能更加復雜，你需要看 機器學習中的數學(4)-線性判別分析（LDA）, 主成分分析(PCA)

反正PCA最終的目的就是讓原本的數據在降維后方差最大，也就是能把數據分開。我們一般人記住這一點，然后會用軟件，能解讀結果就好。

R語言能做主成分分析功能很多，比如說prcomp,princomp，principal. 這里用到是prcomp,它使用的是奇異值分解而不是特征值分解，至于這兩個有啥區(qū)別，其實這不是我關心的，當然你想深入了解到話，建議閱讀 機器學習中的數學(5)-強大的矩陣奇異值分解(SVD)及其應用

代碼如下：

pca.results <- prcomp(iris[1:4],center = TRUE, retx = T)

prcomp會返回一個列表，里面包含如下內容，你可以用print(pca.results)查看

• sdev: 對應每個主成分的標準偏差，值越大說明解釋變異越大。
• rotation: 變量的負荷(loading)矩陣，用于衡量一個變量在各個主成分重要性。
• x: 原矩陣經線性變化后的矩陣（下圖不顯示，可以用pca.results$x提取)

PCA results

雖然有現成的工具提取pca.results里的數據進行作圖，但這樣無助于理解，我們需要自己手動提取結果進行作圖。

iris_rotate_df <- as.data.frame(pca.results$x[,c("PC1","PC2")])
iris_rotate_df$Species <- iris$Species
p <- ggplot(iris_rotate_df, aes(x=PC1, y=PC2)) + geom_point(aes(colour=Species))
print(p)

PCA plot

這個結果雖然也無法完美地區(qū)分出"versicolor"和"virginica"，但是比我們盲目的尋找好多了。這里的主成分是原來的各個變量線性組合結果。

當然我們還有一個問題，各個變量在不同的主成分里的權重是多少？這個就要看負荷矩陣了,也就是Rotation部分。

在PCA中，我們將協方差矩陣分成了標量部分（特征值）和有向部分（特征向量），通過計算特征向量和開方后特征值的乘積，就稱之為負荷(loading)

我們發(fā)現在PC1里面Petal.Length的絕對值最大，而在PC2里面，Sepal.Width的絕對值最大，于是我們就可以從原來的數據集中提取這兩個變量進行作圖了。

ggplot(iris,aes(x=Sepal.Width, y=Petal.Length)) + geom_point(aes(colour=Species))

下面的作圖直接畫圖的結果，是不是感覺和上面的PCA圖很像。當我手動將其旋轉后得到下面的右圖，你會發(fā)現這和PCA的結果幾乎一摸一樣。

loading

那么我們如何在圖上展示各個變量在各個PC的占比呢？這個就需要畫幾個箭頭了，

p + annotate("segment", x=0,xend=0.35,y=0,yend=-0.65,arrow=arrow()) +
  annotate("text", x=0.35,y=-0.68, label="Sepal.Length") +
  annotate("segment", x=0,xend=-0.08,y=0,yend=-0.73,arrow=arrow()) +
  annotate("text", x=-0.08,y=-0.75, label="Sepal.Width") 

下圖中"Septal.Width"朝下，Petal.Length朝右，如果你手動旋轉一下，就是差不多是以"Septal.With"為X軸,"Petal.Length"為Y軸的結果了。 所以箭頭的朝向不重要，重點是長度。

arrow

當然也有現成的包來完成上圖, 不過知道原理后的你應該會解讀結果了吧

library(ggord)
ggord(pca.results, iris$Species, size=1.25)

ggord

當然要想真正的理解PCA分析，你還是了解一點線性代數的知識