主成分分析（PCA：Principal Component Analysis）非常有助於我們理解高維資料，我利用Stack Overflow的每日訪問資料對主成分分析進行了實踐和探索，你可以在rstudio :: conf 2018上找到其中一篇演講的錄音。演講的重點主要是我對於PCA的理解，而這篇文章中，我將主要介紹我是如何實現PCA的，以及我是如何製作演講中使用到的圖表的。
rstudio :: conf 2018
https://www.rstudio.com/resources/videos/understanding-pca-using-shiny-and-stack-overflow-data/
高維資料
此次分析使用的是去年Stack Overflow上註冊使用者訪問量前500的標籤資料。為了簡化處理，本文只使用了10％的註冊流量資料進行分析，但實際上我已經對所有流量資料進行了類似的分析，並獲得了幾乎相同的結果。
標籤資料
https://stackoverflow.com/tags
現在，把每個註冊使用者都想象成高維空間中的一個點，空間的座標軸是R、JavaScript、C++等技術。那麼，在這個高維空間中，做相似工作的人對應的點就會彼此接近。接下來PCA會把這個高維空間轉變成一個新的具有特殊特徵的“特殊”高維空間。
在資料庫中適當地抽取資料後，最開始的資料看起來就像下面這樣：
library(tidyverse) library(scales) tag_percents ## # A tibble: 28,791,663 x 3 ## User Tag Value ## <int> <chr> <dbl> ## 1 1 exception-handling 0.000948 ## 2 1 jsp 0.000948 ## 3 1 merge 0.00284 ## 4 1 casting 0.00569 ## 5 1 io 0.000948 ## 6 1 twitter-bootstrap-3 0.00569 ## 7 1 sorting 0.00474 ## 8 1 mysql 0.000948 ## 9 1 svg 0.000948 ## 10 1 model-view-controller 0.000948 ## # ... with 28,791,653 more rows
可以看出，資料很乾淨，每行只有使用者編號和技術標籤。這裡的User列是隨機ID，而非Stack Overflow的識別符號。在Stack Overflow中，我們公開了大量資料，但流量資料（即哪些使用者訪問過哪些問題）是沒有公開的。
對高維資料進行真正的匿名化其實是非常困難的，而這裡為了進行脫敏處理，我的做法是隨機化資料順序，並用數字替換Stack Overflow的識別符號。Value列表示過去一年該使用者對該標籤的瀏覽量佔該標籤總瀏覽量的比例。
部分資料連結：
https://stackoverflow.blog/2010/06/13/introducing-stack-exchange-data-explorer/
https://cloud.google.com/bigquery/public-data/stackoverflow，
https://meta.stackexchange.com/questions/19579/where-are-the-stack-exchange-data-dumps
先不考慮脫敏的問題，我們首先看看使用者主要瀏覽的技術標籤有哪些，這張圖表給了我們一個直觀的概念。.
tag_percents %>% group_by(Tag) %>% summarise(Value = mean(Value)) %>% arrange(desc(Value)) %>% top_n(15) %>% mutate(Tag = reorder(Tag, Value)) %>% ggplot(aes(Tag, Value, label = Tag, fill = Tag)) + geom_col(alpha = 0.9, show.legend = FALSE) + geom_text(aes(Tag, 0.001), hjust = 0, color = "white", size = 4, family = "IBMPlexSans-Bold") + coord_flip() + labs(x = NULL, y = "Average % of a user's traffic") + scale_y_continuous(labels = percent_format(), expand = c(0.015,0)) + theme(axis.text.y=element_blank())
實施PCA
我們喜歡乾淨的資料，一是因為它就是我們查詢資料庫的結果，二是因為它可用於實現PCA等機器學習演算法的探索性資料分析。為了實現PCA，我們需要一個矩陣，在這個例子裡稀疏矩陣（sparse matrix）就是最佳選擇——因為大多數開發人員只訪問一小部分技術標籤，因此我們的矩陣中會有很多零。tidytext軟體包中有一個函式cast_sparse()，它可以把上面的資料轉換為稀疏矩陣。
sparse_tag_matrix <- tag_percents %>% tidytext::cast_sparse(User, Tag, Value)
R中有幾個實現PCA的演算法是體會不到稀疏矩陣的美感的，比如prcomp()——此演算法的第一步就是將剛剛製作好的稀疏矩陣強制轉換成一個常規矩陣，然後你要在那裡幹坐一輩子等它執行完，因為在它執行的時候電腦根本沒有記憶體讓你去做其他事了（別問我是怎麼知道的）。當然，R中也有一個程式包利用了稀疏矩陣的優勢——irlba。
在建立模型前，也別忘記先用scale()函式將你的矩陣規範化，這對於PCA的實現非常重要。
tags_scaled <- scale(sparse_tag_matrix) tags_pca <- irlba::prcomp_irlba(tags_scaled, n = 64)
其中prcomp_irlba()函式的引數n代表我們想要得到的主成分個數。
那麼這一步究竟發生了什麼？我們會在接下來的章節中慢慢介紹。
class(tags_pca) ## [1] "irlba_prcomp" "prcomp" names(tags_pca) ## [1] "scale" "totalvar" "sdev" "rotation" "center" "x" ?
PCA的結果分析
我喜歡處理資料框格式的資料，所以接下來我要用tidy()函式來整理我的PCA結果，以便用dplyr包處理輸出結果和用ggplot2繪圖。 broom包並不能完美地處理irlba的輸出結果，所以我會將它們與我自己的資料框經過一點修整後合併到一起。
library(broom) tidied_pca <- bind_cols(Tag = colnames(tags_scaled), tidy(tags_pca$rotation)) %>% gather(PC, Contribution, PC1:PC64) tidied_pca ## # A tibble: 39,232 x 3 ## Tag PC Contribution ## <chr> <chr> <dbl> ## 1 exception-handling PC1 -0.0512 ## 2 jsp PC1 0.00767 ## 3 merge PC1 -0.0343 ## 4 casting PC1 -0.0609 ## 5 io PC1 -0.0804 ## 6 twitter-bootstrap-3 PC1 0.0855 ## 7 sorting PC1 -0.0491 ## 8 mysql PC1 0.0444 ## 9 svg PC1 0.0409 ## 10 model-view-controller PC1 0.0398 ## # ... with 39,222 more rows
注意到這裡我的資料框的每一行只有一個技術標籤及它構成的主成分。
那麼從整體來看，這些結果又是什麼樣子的呢？請見下圖：
tidied_pca %>% filter(PC %in% paste0("PC", 1:6)) %>% ggplot(aes(Tag, Contribution, fill = Tag)) + geom_col(show.legend = FALSE, alpha = 0.8) + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + labs(x = "Stack Overflow tags", y = "Relative importance in each principal component") + facet_wrap(~ PC, ncol = 2)
很漂亮吧有木有！我們上面看的是前六個主成分，圖中x軸上是按字母順序排列的單個Stack Overflow標籤，縱軸表示該技術標籤對這一PC的貢獻度。我們也可以看出有關聯的技術可能是以相同的字母開頭，故而會排列在一起，例如PC4中的橙色等。
下面讓我們主要分析一下第一個主成分的構成。
tidied_pca %>% filter(PC == "PC1") %>% top_n(40, abs(Contribution)) %>% mutate(Tag = reorder(Tag, Contribution)) %>% ggplot(aes(Tag, Contribution, fill = Tag)) + geom_col(show.legend = FALSE, alpha = 0.8) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), xis.ticks.x = element_blank()) + labs(x = "Stack Overflow tags", y = "Relative importance in principle component")
現在我們可以看到哪些技術標籤對這個成分有貢獻。從貢獻為正的標籤來看，主要有前端Web開發技術，如HTML、JavaScript、jQuery、CSS等。從貢獻為負的標籤來看，主要有Python，C ++以及低階技術詞彙，如字串（strings）、列表（lists）等。這意味著Stack Overflow的使用者之間最大的差異在於他們是使用前端Web技術更多一些還是Python和一些低階技術更多一些。
那麼第二個主成分又是怎樣的呢？
tidied_pca %>% filter(PC == "PC2") %>% top_n(40, abs(Contribution)) %>% mutate(Tag = reorder(Tag, Contribution)) %>% ggplot(aes(Tag, Contribution, fill = Tag)) + geom_col(show.legend = FALSE, alpha = 0.8) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.ticks.x = element_blank()) + labs(x = "Stack Overflow tags", y = "Relative importance in principle component")
第一個主成分是兩種軟體工程的對比，但第二個主成分則更像是一個結果為是/否的二分類變數。它告訴了我們開發人員工作中是否使用C＃、.NET、Visual Studio和Microsoft技術堆疊的其餘部分。這意味著Stack Overflow的使用者之間的第二大差異在於他們是否訪問了這些型別的微軟技術問題。
我們可以繼續研究其他的主成分，瞭解更多關於Stack Overflow技術生態系統的知識，但其實我已經在視訊中進行了相關內容的講解，也研究了那些與我們資料科學人員相關的技術。我還製作了一個名叫Shiny的應用程式，在上面你可以隨意選擇你想研究的主成分。而且我敢打賭，只要你用過一次Shiny，你就能想象到我是如何開始這項研究的！
高維平面的對映
PCA最酷的地方在於它能幫我們思考和推理高維資料，其中一項功能就是將高維資料對映到可繪圖的二維平面上。接下來我們來看看它是如何做到這一點的。
其實這一步用broom :: augment()就能實現，並且還能計算出每個成分對整個資料集方差解釋的百分比。
percent_variation <- tags_pca$sdev^2 / sum(tags_pca$sdev^2) augmented_pca <- bind_cols(User = rownames(tags_scaled), tidy(tags_pca$x)) augmented_pca ## # A tibble: 164,915 x 65 ## User PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 ## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 1 2.16 5.70 1.63 0.967 0.0214 -1.37 -1.98 -2.94 -0.860 ## 2 2 0.350 3.38 -6.12 -10.0 1.39 0.882 5.35 -3.30 -2.73 ## 3 3 2.75 -3.91 0.801 1.73 1.24 -0.837 2.03 2.76 0.300 ## 4 4 3.27 -3.37 -1.00 2.39 -3.59 -2.68 0.449 -2.82 -1.25 ## 5 5 9.44 -4.24 3.88 -1.62 -2.96 4.01 -1.32 -3.54 3.25 ## 6 6 5.47 -5.13 1.57 2.94 -0.170 0.342 3.34 6.09 1.72 ## 7 7 4.30 -0.442 -1.52 0.329 -2.13 0.908 -3.30 -5.02 -1.39 ## 8 8 -0.691 0.668 -1.76 -7.74 -2.94 -5.28 -9.71 5.28 0.732 ## 9 9 3.84 -2.65 0.760 1.34 2.06 -0.927 1.35 5.11 -2.69 ## 10 10 3.23 4.13 2.81 2.68 -1.12 -1.30 -0.319 -1.23 -0.723 ## # ... with 164,905 more rows, and 55 more variables: PC10 <dbl>, ## # PC11 <dbl>, PC12 <dbl>, PC13 <dbl>, PC14 <dbl>, PC15 <dbl>, ## # PC16 <dbl>, PC17 <dbl>, PC18 <dbl>, PC19 <dbl>, PC20 <dbl>, ## # PC21 <dbl>, PC22 <dbl>, PC23 <dbl>, PC24 <dbl>, PC25 <dbl>, ## # PC26 <dbl>, PC27 <dbl>, PC28 <dbl>, PC29 <dbl>, PC30 <dbl>, ## # PC31 <dbl>, PC32 <dbl>, PC33 <dbl>, PC34 <dbl>, PC35 <dbl>, ## # PC36 <dbl>, PC37 <dbl>, PC38 <dbl>, PC39 <dbl>, PC40 <dbl>, ## # PC41 <dbl>, PC42 <dbl>, PC43 <dbl>, PC44 <dbl>, PC45 <dbl>, ## # PC46 <dbl>, PC47 <dbl>, PC48 <dbl>, PC49 <dbl>, PC50 <dbl>, ## # PC51 <dbl>, PC52 <dbl>, PC53 <dbl>, PC54 <dbl>, PC55 <dbl>, ## # PC56 <dbl>, PC57 <dbl>, PC58 <dbl>, PC59 <dbl>, PC60 <dbl>, ## # PC61 <dbl>, PC62 <dbl>, PC63 <dbl>, PC64 <dbl>
注意到這裡我其實有更廣闊的資料框可供使用，並且我還沒有使用gather()函式——為了便於繪圖。物件percent_variation是一個向量，它包含了每個主成分對整個資料集的方差解釋的百分比。
augmented_pca %>% mutate(User = as.integer(User)) %>% filter(User %% 2 == 0) %>% ggplot(aes(PC1, PC2)) + geom_point(size = 1.3, color = "midnightblue", alpha = 0.1) + labs(x = paste0("Principal component 1 (", percent(percent_variation[1]), ")"), y = paste0("Principal component 2 (", percent(percent_variation[2]),")"), title = "Projection of Stack Overflow traffic on to the first two principal components", subtitle = "The very high dimensional space can be projected down onto components we have explored")
可以看出，為了儘量減少過度繪圖，這個圖裡我把每兩個人用一個點表示。還記得第一個主成分是前端開發人員到Python和低階技術人員的橫向擴充，而第二個主成分則全部是關於微軟技術堆疊的。由上我們可以看到描述Stack Overflow標籤的高維資料是如何投影到前兩個主成分的。可以注意到我已在每個軸中新增了方差百分比，同時這些數字並不是很高，這也與我們現實生活中的情況相吻合，即事實上Stack Overflow的使用者之間差異很大，如果你想將這些主成分中的任意一個用於降維或作為模型中的預測變數，請慎重考慮。
應用
說到現實生活，我發現PCA非常有助於我們理解高維資料集。比如說，基於完全相同的資料，我最近在使用PCA探索的另一個問題是亞馬遜可能考慮讓哪些城市成為其第二總部。實際上，PCA給出的主成分結果以及不同技術對其的貢獻率已經不盡相同——因為幾個月已經過去了，而且使用者們在高維空間中也不是完全靜止的。如果你有任何問題或反饋，請及時聯絡我。
相關報導：
https://juliasilge.com/blog/stack-overflow-pca/