手把手 | 用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
##    <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/

原文发布于微信公众号 - 大数据文摘(BigDataDigest)

原文发表时间:2018-05-31

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏我和未来有约会

[Silverlight动画]转向行为 - 群落

说到群落,很难不引用Craig Reynolds和他的"boilds"模拟系统。Reynolds很牛的将一个看似非常恐怖的复杂过程,拆成了几个比较简单的行为。 ...

2258
来自专栏WOLFRAM

用 Wolfram 语言制作圣诞动画

2212
来自专栏PPV课数据科学社区

如何向小白介绍何谓机器学习和数据挖掘?买回芒果他就懂了

  买芒果   嘴馋的你想吃芒果了,于是你走到水果摊,挑了几个让老板过过秤,然后你再根据芒果的斤两付钱走人。   显然,买芒果你当然是挑着最甜、最熟的来买(因为...

3607
来自专栏开发 & 算法杂谈

动态数据竞争检测方法实验分析(一)

之前的文章大致介绍了一下我们的动态数据竞争检测平台如何构建,这篇文章主要是在动态数据竞争检测平台上实现了之前介绍的数据竞争检测方法,我们扩展了其中的一些方法使得...

2272
来自专栏余林丰

初学数据挖掘——相似性度量(二)

  上一篇中介绍了四个算法,并用四个算法分别计算了两个人的相似度。这篇就来讲讲相似性算法在实际当中怎么用。第一:将指定的人与其他人作相似性比较,并从高到低进行排...

2436
来自专栏新智元

邓侃:谷歌Talk to books引爆搜索方式革命

?---- 【新智元导读】昨天,新智元介绍了谷歌的全新搜索工具“Talk to Books”,基于自然语言文本理解,用户能够凭语义而非关键词来实现搜索功能。谷歌...

2966
来自专栏窗户

怀念Galois

  我的第一篇谈到具体学科的博客,还是献给我最钟爱的数学。   个人比较喜欢离散数学,并非因为曲高和寡,而是因为数学分析、概率论、拓扑学、泛函之类的高手实在太多...

1965
来自专栏数据科学与人工智能

【机器学习】机器学习大白话

买芒果 嘴馋的你想吃芒果了,于是你走到水果摊,挑了几个让老板过过秤,然后你再根据芒果的斤两付钱走人。 显然,买芒果你当然是挑着最甜、最熟的来买(因为你是根据重量...

2615
来自专栏奇点大数据

遗传算法(2)

在遗传算法中我们再举一个求极大值的例子。这种例子也是比较多见的,只要我们把一些数据关系描述成函数之后就会有一些求极大值或者极小值的问题。 其实极大值和极小值是一...

35012
来自专栏hadoop学习笔记

Hanlp自然语言处理工具的使用演练

Hanlp是由一系列模型与算法组成的工具包,目标是普及自然语言处理在生产环境中的应用。Hanlp具备功能完善、性能高效、架构清洗、语料时新、可自定义的特点;提供...

2032

扫码关注云+社区

领取腾讯云代金券