05.31 用StackOverflow訪問數據實現主成分分析（PCA）

大數據文摘出品

編譯：汪小七、張馨月、雲舟

主成分分析（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
##   
## 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 %
tidytext::cast_sparse(User, Tag, Value)

R中有幾個實現PCA的算法是體會不到稀疏矩陣的美感的，比如prcomp()——此算法的第一步就是將剛剛製作好的稀疏矩陣強制轉換成一個常規矩陣，然後你要在那裡乾坐一輩子等它運行完，因為在它運行的時候電腦根本沒有內存讓你去做其他事了（別問我是怎麼知道的）。當然，R中也有一個程序包利用了稀疏矩陣的優勢——irlba。

在建立模型前，也別忘記先用scale()函數將你的矩陣規範化，這對於PCA的實現非常重要。

tags_scaled tags_pca 

其中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  tidy(tags_pca$rotation)) %>%
 gather(PC, Contribution, PC1:PC64)
tidied_pca
## # A tibble: 39,232 x 3
## Tag PC Contribution
##   
## 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 augmented_pca  tidy(tags_pca$x)) 

augmented_pca
## # A tibble: 164,915 x 65
## User PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
##          
## 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 ,
## # PC11 , PC12 , PC13 , PC14 , PC15 ,
## # PC16 , PC17 , PC18 , PC19 , PC20 ,
## # PC21 , PC22 , PC23 , PC24 , PC25 ,
## # PC26 , PC27 , PC28 , PC29 , PC30 , 

## # PC31 , PC32 , PC33 , PC34 , PC35 ,
## # PC36 , PC37 , PC38 , PC39 , PC40 ,
## # PC41 , PC42 , PC43 , PC44 , PC45 ,
## # PC46 , PC47 , PC48 , PC49 , PC50 ,
## # PC51 , PC52 , PC53 , PC54 , PC55 ,
## # PC56 , PC57 , PC58 , PC59 , PC60 ,
## # PC61 , PC62 , PC63 , PC64 
 
 

注意到這裡我其實有更廣闊的數據框可供使用，並且我還沒有使用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/