[原] 数据科学教程：R语言与DataFrame[2019版]

前言

9102年是互联网大环境不太好的一年，这一年更需要苦练基本功，数据科学领域的基本功无非就是数据处理，而 DataFrame 是其中的核心。那么，都9102年了，如何在 R 语言中优雅地使用 DataFrame 呢？是否优雅是矿工生产力差异的重要来源，本文将介绍最近三年来 DataFrame 的最近进展。

R 中的 DataFrame 新趋势

1. tidy 编程流
2. 分布式计算
3. 集成 NoSQL
4. 集成各领域核心算法

按照功能来划分：

1. 核心计算：tidyverse/sparklyr
2. 图计算： tidygraph/graphframes
3. 空间计算：sf/geospark/tidyraster
4. 时间序列：tidyquant
5. 文本计算：tidytext/tabulizer
6. 模型构建：tidymodels
7. 仿真建模：simmer
8. 贝叶斯分析：tidybayes
9. 深度学习：keras

核心计算

tidyverse 介绍

谈到数据处理的基本功必须先了解 DataFrame 的最核心能力：比SQL还要优雅地处理数据。作为第四代数据处理语言的核心 tidyverse 是必须要了解的。

tidyverse 实现了核心计算的大一统，它囊括了八个数据科学最高频的操作：

• ggplot2, 支持数据可视化
• dplyr, 数据操作
• tidyr, 数据清洗
• readr, 数据IO
• purrr, 函数式编程
• tibble, 数据框轻量化格式
• stringr, 文本处理
• forcats, 因子类型处理

对于 tidyverse 的各个子模块可以参考 R语言与DataFrame 2016版, 其中详细介绍了各模块的情况，本文不再赘述，更多着墨与新趋势。八大模块中，dplyr 是核心中的核心，因为它实现了数据处理的前后端分离，支持 R 语言的前端和多种结构化数据库的后端，比如 RPostGIS,RMySQL,sparklyr 等。

Spark 数据处理举例

下面举一个 sparklyr 处理 iris 数据的例子来感受一下现代数据处理：

步骤一：数据库连接

library(sparklyr)
library(tidyverse)

# 分布式OLAP数据库 Hive 连接
sc <- sparklyr::spark_connect(master = "localhost",
                             spark_home = "/data/FinanceR/Spark",
                             version = "2.2.0",
                             config = sparklyr::spark_config())

# 或者分布式OLTP数据库 tidb 连接
# sc <- RMySQL::dbConnect(RMySQL::MySQL(),
#         user="FinanceR", password="FinanceR",
#         dbname="FinanceR", host="192.168.1.100")

# 或者基于内存索引的 MPP数据库 impala 链接
# impala <- implyr::src_impala(
#    drv = odbc::odbc(),
#    driver = "Cloudera ODBC Driver for Impala",
#    host = "host",
#    port = 21050,
#    database = "default",
#    uid = "username",
#    pwd = "password"
#  )

步骤二：数据操作

# 在 mutate 中支持 Hive UDF

remote_df = copy_to(sc,iris,"iris_tbl")

remote_df = dplyr::tbl(sc,from = "db.iris_tbl") # 定义数据源表 

# 或者 remote_df = dplyr::tbl(sc,from = dplyr::sql("select * from db.iris_tbl limit 10")) #


# 以比SQL更接近于自然语言的方法编程，以基本的 过滤、聚合、排序为例

remote_df %>%
    dplyr::mutate(a = Sepal_Length+2) %>%   # 在 mutate 中支持 Hive UDF
    dplyr::filter(a > 2)%>%
    dplyr::group_by(Species)%>%
    dplyr::summarize(count = n())%>%
    dplyr::select(cnt = count)%>% 
    dplyr::arrange(desc(cnt)) %>%
    dplyr::mutate_if(is.double,as.character) -> # 配合 逆天的批量操作
    pipeline

步骤三：检验语法

## 一顿操作猛如虎之后，将数据操作的对应SQL执行过程打印出来

pipeline %>% 
    dbplyr::sql_render() 
    
# <SQL> SELECT CAST(`cnt` AS STRING) AS `cnt`
FROM (SELECT *
FROM (SELECT `count` AS `cnt`
FROM (SELECT `Species`, count(*) AS `count`
FROM (SELECT *
FROM (SELECT `Sepal_Length`, `Sepal_Width`, `Petal_Length`, `Petal_Width`, `Species`, `Sepal_Length` + 2.0 AS `a`
FROM `iris_tbl`) `nmpfznbuzf`
WHERE (`a` > 2.0)) `ygjaktrzzu`
GROUP BY `Species`) `lbmpsjksox`) `tdrwlnidxw`
ORDER BY `cnt` DESC) `mzuwuocxay`

步骤四：数据存储

## 将 pipeline 结果自动缓存为临时表 tmp_tbl    
tmp_tbl = pipeline %>% compute("tmp_tbl")

# 将 临时表按照条件过滤后再写回 Hive 数据库，省去建表过程

tmp_tbl %>% 
   filter(cnt > 10) %>%
   sparklyr::spark_write_table("db.financer_res")

通过上面的例子，我们可以看到使用 sparklyr 而不是直接使用 Hive ，能够使得数据分析流程变得更加高效和方便。

空间计算

在 Edzer 等人的努力之下，空间计算社区逐渐建立起来一套优雅的以 sf 为核心的空间计算工作流。

核心组成主要包括

1. sf：矢量空间核心计算框架，实现九交模型、Voronoi图、凸包计算、GROUP BY、数据IO等标准操作。
2. stars：栅格空间核心计算框架，实现了核心multi-dimension的数据结构和数据IO，使得卫星影像时序数据处理更加高效。
3. rpostgisLT：基于PostGIS的轨迹数据分析工具
4. leaflet：轻量级的空间可视化组件
5. mapview：基于leaflet的高级空间可视化集成方案
6. mapedit：基于leaflet和shiny的空间元素编辑工具
7. tmap：专题图的绘制
8. stplanr：可持续交通规划

在空间数据处理中，通常将数据分为两大类：

1. 矢量空间数据（连续）
2. 栅格空间数据（离散）

下面按照空间数据类型分别介绍：

矢量空间数据处理：sf

sf 是 OGC 和 ISO 共同拟定的一套通用二维几何图形(包括点, 线, 面, 多点, 多线 等）存储和访问GIS模型, 旨在标准化要素存取流程。sf 的简称是 Simple Feature Access, 意为 简单要素存取。

空间数据格式转化

在 spatial tidy workflow 中的核心就是 sf 的空间数据结构。

以 trip network 举例，实现 data.frame 到 sf 的空间 dataframe 的转化：

library(tidyverse)
library(stplanr)
data(cents, flow)
travel_network <- stplanr::od2line(flow = flow, zones = cents) %>%  
                  sf::st_as_sf() %>%
                  sf::st_set_transform(4326)

travel_network %>% 
    select( Area.of.residence)

Simple feature collection with 49 features and 1 field
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -1.550806 ymin: 53.8041 xmax: -1.511861 ymax: 53.82887
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
       Area.of.residence                       geometry
920573         E02002361 LINESTRING (-1.516734 53.82...
920575         E02002361 LINESTRING (-1.516734 53.82...
920578         E02002361 LINESTRING (-1.516734 53.82...
920582         E02002361 LINESTRING (-1.516734 53.82...
920587         E02002361 LINESTRING (-1.516734 53.82...
920591         E02002361 LINESTRING (-1.516734 53.82...
920601         E02002361 LINESTRING (-1.516734 53.82...
921220         E02002363 LINESTRING (-1.535617 53.82...
921222         E02002363 LINESTRING (-1.535617 53.82...
921225         E02002363 LINESTRING (-1.535617 53.82...

从 travel_network 对象中可以看到，sf 实现了 空间信息（坐标系、集合类型、维度、bbox、投影系） + DataFrame 的完美结合, 支持与 dplyr 操作的完美结合，使得空间数据的处理变得简单优雅。

以骑行网络上自行车出行的流量分布统计为例，通过 sf::st_set_geometry(NULL) 函数将空间信息隐去：

travel_network %>%
 dplyr::group_by(Bicycle) %>% 
 dplyr::summarise(cnt = n()) %>% 
 sf::st_set_geometry(NULL)

# A tibble: 9 x 2
  Bicycle   cnt
*   <int> <int>
1       0    24
2       1    11
3       2     8
4       3     1
5       4     1
6       5     1
7       6     1
8      10     1
9      12     1

空间数据可视化

travel_network 对象通过 leaflet 直接实现快速可视化，不需要额外指定经纬度参数：

travel_network %>%
leaflet::leaflet() %>%
leaflet::addTiles() %>%
leaflet::addPolylines()

sf 的架构分析

sf 本质上是若干 C++/Rcpp库的一个封装，旨在为空间数据分析提供一个统一的工作流接口，下面是 sf 的框架图，包括以下几个核心构成：

1. tidyverse: 核心数据操作包括 join/select/group_by/filter/arrange 等
2. rgdal: 多种格式的数据IO
3. rgeos: 基本空间计算算子，主要包括二维空间操作
4. rproj4：投影系、坐标系的选择与变换
5. lwgeom：高级空间计算算子，主要包括球面空间操作、PostGIS 插件等
6. units：计量单位处理，比如 米、度 等
7. DBI：关系型数据库交互的核心接口比如 RPostGIS

sp 与 sf 对比操作

sf 的核心数据操作主要是基于 sp 包实现的，sp 包的api由于历史原因实现的功能非常丰富但是命名并不规则，sf 在它的基础上对接了 Simple Feature Access 标准。下面是 sp 操作 和 sf 操作的对比

辅助操作

数据格式转化

空间几何计算

九交模型操作

维度扩展九交模型，旨在解决二维空间的集合拓扑关系判断。维度扩展九交模型的形式如下所示，假设 A、B 表示两个相互独立的多边形，则有 I（A）和I（B）表示A和B的内部（inside），B（A）和B（B）表示A和B的边界（border），E（A)和E（B）表示A和B的外部。下面是DE-9IM 模型的示意图：

通过九交模型可以支持完整的二维空间八类拓扑关系，对应的 SF 函数如下：

高级操作

更多空间计算操作 详见 SF CheatSheet.

GeoSpark 数据处理举例

GeoSpark 解决了大规模矢量数据的空间索引问题，进而产生了比如大规模处理 九交模型 的对应方案，解决了PG/SF单机处理矢量空间数据的性能瓶颈问题。

以 DISJOINT 为例：

步骤一：环境初始化

连接 Spark 程序

pak::pkg_install("harryprince/geospark")

library(sparklyr)
library(dplyr)
library(geospark)
conf = spark_config()
conf$`sparklyr.instances.num` = 50
conf$spark.serializer <- "org.apache.spark.serializer.KryoSerializer"
conf$spark.kryo.registrator <- "org.datasyslab.geospark.serde.GeoSparkKryoRegistrator"
Sys.set(SPARK_HOME=”~/spark”)
sc <- spark_connect(master = "yarn-client",config = conf)

注册空间算子

register_gis(sc)

步骤二：

本地伪造数据

polygons <- read.table(text="california area|POLYGON ((-126.4746 32.99024, -126.4746 42.55308, -115.4004 42.55308, -115.4004 32.99024, -126.4746 32.99024))
new york area|POLYGON ((-80.50781 36.24427, -80.50781 41.96766, -70.75195 41.96766, -70.75195 36.24427, -80.50781 36.24427))
texas area |POLYGON ((-106.5234 25.40358, -106.5234 36.66842, -91.14258 36.66842, -91.14258 25.40358, -106.5234 25.40358))
dakota area|POLYGON ((-106.084 44.21371, -106.084 49.66763, -95.71289 49.66763, -95.71289 44.21371, -106.084 44.21371))
", sep="|",col.names=c("area","geom"))

points <- read.table(text="New York|NY|POINT (-73.97759 40.74618)
New York|NY|POINT (-73.97231 40.75216)
New York|NY|POINT (-73.99337 40.7551)
West Nyack|NY|POINT (-74.06083 41.16094)
West Point|NY|POINT (-73.9788 41.37611)
West Point|NY|POINT (-74.3547 41.38782)
Westtown|NY|POINT (-74.54593 41.33403)
Floral Park|NY|POINT (-73.70475 40.7232)
Floral Park|NY|POINT (-73.60177 40.75476)
Elmira|NY|POINT (-76.79217 42.09192)
Elmira|NY|POINT (-76.75089 42.14728)
Elmira|NY|POINT (-76.84497 42.12927)
Elmira|NY|POINT (-76.80393 42.07202)
Elmira|NY|POINT (-76.83686 42.08782)
Elmira|NY|POINT (-76.75089 42.14728)
Alcester|SD|POINT (-96.63848 42.97422)
Aurora|SD|POINT (-96.67784 44.28706)
Baltic|SD|POINT (-96.74702 43.72627)
Beresford|SD|POINT (-96.79091 43.06999)
Brandon|SD|POINT (-96.58362 43.59001)
Minot|ND|POINT (-101.2744 48.22642)
Abercrombie|ND|POINT (-96.73165 46.44846)
Absaraka|ND|POINT (-97.21459 46.85969)
Amenia|ND|POINT (-97.25029 47.02829)
Argusville|ND|POINT (-96.95043 47.0571)
Arthur|ND|POINT (-97.2147 47.10167)
Ayr|ND|POINT (-97.45571 47.02031)
Barney|ND|POINT (-96.99819 46.30418)
Blanchard|ND|POINT (-97.25077 47.3312)
Buffalo|ND|POINT (-97.54484 46.92017)
Austin|TX|POINT (-97.77126 30.32637)
Austin|TX|POINT (-97.77126 30.32637)
Addison|TX|POINT (-96.83751 32.96129)
Allen|TX|POINT (-96.62447 33.09285)
Carrollton|TX|POINT (-96.89163 32.96037)
Carrollton|TX|POINT (-96.89773 33.00542)
Carrollton|TX|POINT (-97.11628 33.20743)
Celina|TX|POINT (-96.76129 33.32793)
Carrollton|TX|POINT (-96.89328 33.03056)
Carrollton|TX|POINT (-96.77763 32.76727)
Los Angeles|CA|POINT (-118.2488 33.97291)
Los Angeles|CA|POINT (-118.2485 33.94832)
Los Angeles|CA|POINT (-118.276 33.96271)
Los Angeles|CA|POINT (-118.3076 34.07711)
Los Angeles|CA|POINT (-118.3085 34.05891)
Los Angeles|CA|POINT (-118.2943 34.04835)
Los Angeles|CA|POINT (-118.2829 34.02645)
Los Angeles|CA|POINT (-118.3371 34.00975)
Los Angeles|CA|POINT (-118.2987 33.78659)
Los Angeles|CA|POINT (-118.3148 34.06271)", sep="|",col.names=c("city","state","geom"))

points <- read.table(text="New York|NY|POINT (-73.97759 40.74618)
New York|NY|POINT (-73.97231 40.75216)
New York|NY|POINT (-73.99337 40.7551)
West Nyack|NY|POINT (-74.06083 41.16094)
West Point|NY|POINT (-73.9788 41.37611)
West Point|NY|POINT (-74.3547 41.38782)
Westtown|NY|POINT (-74.54593 41.33403)
Floral Park|NY|POINT (-73.70475 40.7232)
Floral Park|NY|POINT (-73.60177 40.75476)
Elmira|NY|POINT (-76.79217 42.09192)
Elmira|NY|POINT (-76.75089 42.14728)
Elmira|NY|POINT (-76.84497 42.12927)
Elmira|NY|POINT (-76.80393 42.07202)
Elmira|NY|POINT (-76.83686 42.08782)
Elmira|NY|POINT (-76.75089 42.14728)
Alcester|SD|POINT (-96.63848 42.97422)
Aurora|SD|POINT (-96.67784 44.28706)
Baltic|SD|POINT (-96.74702 43.72627)
Beresford|SD|POINT (-96.79091 43.06999)
Brandon|SD|POINT (-96.58362 43.59001)
Minot|ND|POINT (-101.2744 48.22642)
Abercrombie|ND|POINT (-96.73165 46.44846)
Absaraka|ND|POINT (-97.21459 46.85969)
Amenia|ND|POINT (-97.25029 47.02829)
Argusville|ND|POINT (-96.95043 47.0571)
Arthur|ND|POINT (-97.2147 47.10167)
Ayr|ND|POINT (-97.45571 47.02031)
Barney|ND|POINT (-96.99819 46.30418)
Blanchard|ND|POINT (-97.25077 47.3312)
Buffalo|ND|POINT (-97.54484 46.92017)
Austin|TX|POINT (-97.77126 30.32637)
Austin|TX|POINT (-97.77126 30.32637)
Addison|TX|POINT (-96.83751 32.96129)
Allen|TX|POINT (-96.62447 33.09285)
Carrollton|TX|POINT (-96.89163 32.96037)
Carrollton|TX|POINT (-96.89773 33.00542)
Carrollton|TX|POINT (-97.11628 33.20743)
Celina|TX|POINT (-96.76129 33.32793)
Carrollton|TX|POINT (-96.89328 33.03056)
Carrollton|TX|POINT (-96.77763 32.76727)
Los Angeles|CA|POINT (-118.2488 33.97291)
Los Angeles|CA|POINT (-118.2485 33.94832)
Los Angeles|CA|POINT (-118.276 33.96271)
Los Angeles|CA|POINT (-118.3076 34.07711)
Los Angeles|CA|POINT (-118.3085 34.05891)
Los Angeles|CA|POINT (-118.2943 34.04835)
Los Angeles|CA|POINT (-118.2829 34.02645)
Los Angeles|CA|POINT (-118.3371 34.00975)
Los Angeles|CA|POINT (-118.2987 33.78659)
Los Angeles|CA|POINT (-118.3148 34.06271)", sep="|",col.names=c("city","state","geom"))

通过 sf 和 mapview 进行 dataframe 快速可视化

M1 = polygons %>% 
sf::st_as_sf(wkt="geom") %>%
 mapview::mapview()


M2 = points %>%
sf::st_as_sf(wkt="geom") %>%
 mapview::mapview()

M1+M2

步骤三：geospark 核心计算

将数据复制到 Spark 集群

polygons_tbl <- copy_to(sc,polygons)
points_tbl <- copy_to(sc, points)

通过 Spark GIS SQL 计算

ex2 <- copy_to(sc,tbl(sc, sql("  SELECT  area,state,count(*) cnt from                            (select area,ST_GeomFromWKT(polygons.geom ,'4326') as y  from polygons)  polygons,                       (SELECT ST_GeomFromWKT (points.geom,'4326') as x,state,city from points) points          where  ST_Contains(polygons.y,points.x) group by area,state")),"test2")

Res = collect(ex2)

关联计算结果，通过 leaflet 展示

Idx_df = polygons %>% 
left_join(Res,by = (c("area"="area"))) %>% 
sf::st_as_sf(wkt="geom")

Idx_df %>% 
leaflet::leaflet() %>% 
leaflet::addTiles() %>% 
leaflet::addPolygons(popup = ~as.character(cnt),color=~colormap::colormap_pal()(cnt))

栅格空间数据处理：stars

stars 建立在 sf 、raster、 tidyverse 之上，旨在简化 Multidimension 数据处理的流程。以公里网格的交通流量统计为例，实现 space x space x mode x time x time x metrics 的多维数据结构：

步骤一： 读取空间数据

library(stars)
## Loading required package: abind
## Loading required package: sf
## Linking to GEOS 3.5.0, GDAL 2.2.2, PROJ 4.8.0
library(dplyr)

nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))

步骤二： 多维数据组合

to = from = st_geometry(nc) # 100 polygons: O and D regions
mode = c("car", "bike", "foot") # travel mode
day = 1:100 # arbitrary
library(units)
units(day) = make_unit("days since 2015-01-01")
hour = set_units(0:23, h) # hour of day
dims = st_dimensions(origin = from,
                         destination = to, 
          mode = mode,
           day = day,
           hour = hour)

图计算

图数据是一种比较特殊的数据结构，通常用 G(V,E) 来表示一个图结构，它由节点(Vertex)和边(Edge)两部分组成。

tidygraph 介绍

tidygraph 提供了清爽的 graph/network 数据操作接口，由 Thomas Lin Pedersen 大神开发，终结了 igraph 和 networkx 的混沌。它使得图结构可以在节点表-边表模式和图模式之间相互切换，一方面，集成了主流图模型，比如 pagerank，louvain 等。另一方面，支持主流图格式，包括 network，phylo，dendrogram，data.tree，graph ，简直是知识图谱和风控领域居家旅行必备神器。

tidygraph 中核心操作见如下：

• 表/图格式转化：as_tibble/as_tbl_graph
• 开启操作节点或边属性操作模式：activate(nodes)/activate(edges)
• 图扩展：bind_edges(),bind_nodes(),bind_graphs()
• 图关联：graph_join()
• 图结构重路由：reroute()
• 图变形/反变形：morph()/ unmorph()
• 图可视化：ggraph
• 中心度指标：centrality_*
• 局部指标：local_*
• 图算子：group_*
• 节点算子：node_*
• 边算子: edge_*
• 遍历算子：bfs_, dfs_, map_local_, map_bfs__，map_dfs__*

以前文提到的 travel_network 为例，将空间结构转化为图结构，并通过 PageRank 计算节点中心度：

library(tidygraph)
library(ggraph)
travel_graph <- sfnet::sfnetworks(travel_network) %>%
                tidygraph::activate(nodes) %>% 
                mutate(degree = tidygraph::centrality_pagerank())

travel_graph

可以看到 travel_graph 生成了以 Node DataFrame + Edge DataFrame 的组合形式，分别记录。
在 active(nodes) 的条件下，优先展示 Node 相关属性。

通过 ggraph 进行可视化呈现

ggraph(travel_graph, layout = 'kk') + 
    geom_edge_fan(aes(alpha = ..index..), show.legend = FALSE) + 
    geom_node_point(aes(size = degree))

支持的算法列表：

中心度指标

• centrality_alpha()
• centrality_authority()
• centrality_betweenness()
• centrality_power()
• centrality_closeness()
• centrality_eigen()
• centrality_hub()
• centrality_pagerank()
• centrality_subgraph()
• centrality_degree()
• centrality_edge_betweenness()

局部指标

• local_size(order = 1, mindist = 0)
• local_members(order = 1, mindist = 0)
• local_triangles()
• local_ave_degree()
• local_transitivity()

图模型

• group_components()
• group_edge_betweenness()
• group_fast_greedy()
• group_infomap()
• group_label_prop()
• group_leading_eigen()
• group_louvain()
• group_optimal()
• group_spinglass()
• group_walktrap()
• group_biconnected_component()

GraphFrames 分布式图模型举例

在大型复杂网络中，GraphFrames 是 tidygraph 很好的替代解决方案，依然以 PageRank 为例：

连接 Spark

library(graphframes)
library(sparklyr)
library(dplyr)

sc <- spark_connect(master = "local", version = "2.3.0")

获取图结构用例

g <- gf_friends(sc)

计算 PageRank

results <- gf_pagerank(g, tol = 0.01, reset_probability = 0.15)
results

GraphFrame 的计算结果：

## GraphFrame
## Vertices:
##   $ id       <chr> "f", "b", "g", "a", "d", "c", "e"
##   $ name     <chr> "Fanny", "Bob", "Gabby", "Alice", "David", "Charlie",...
##   $ age      <int> 36, 36, 60, 34, 29, 30, 32
##   $ pagerank <dbl> 0.3283607, 2.6555078, 0.1799821, 0.4491063, 0.3283607...
## Edges:
##   $ src          <chr> "b", "c", "d", "e", "a", "a", "e", "f"
##   $ dst          <chr> "c", "b", "a", "f", "e", "b", "d", "c"
##   $ relationship <chr> "follow", "follow", "friend", "follow", "friend",...
##   $ weight       <dbl> 1.0, 1.0, 1.0, 0.5, 0.5, 0.5, 0.5, 1.0

总结

过去的三年，是数据科学突飞猛进的三年，大数据和人工智能经历了资本的高潮迭起之后必将回归本质，苦练数据科学基本功，精通数据集的高效清洗与理解才是决胜数据科学的王道。

参考资料

• 数据科学教程：R语言与DataFrame 2016版
• 深入对比数据科学工具箱：Python和R 非结构化数据的结构化
• tidygraph workshop
• https://statkclee.github.io/n...
• https://www.jessesadler.com/p...
• The tidy tools manifesto
• RStudio Spark/Leaflet 与 GIS 最佳实践
• https://www.r-consortium.org/...
• 九交模型原理
• awesome-rmarkdown
• awesome-sparklyr