2 数据获取

数据获取包含两层意思，其一是数据收集，其二是数据搜集。数据收集，往往意味着自己做实验设计，执行实验，回收数据，掌握第一手资料。而数据搜集，往往意味着自己从各个地方搜罗数据，再清洗整理校验，得到可靠的二手或三手数据。从前，统计学家下到试验田，在不同NPK（氮肥、磷肥和钾肥）配比的情况下，收集小麦的产量数据，以确定最佳配比。如今，许多互联网公司都有自己的 App，不断调整产品的结构、排列和色彩，通过 App 收集用户及其行为数据，再以一定的数据模型加工整理成可用的二维表格，分析数据后确定最佳的用户体验。此外，许多政府和非政府的组织机构网站也发布大量的数据，比如各个国家的国家统计局和地方统计局，世界银行，国际货币基金组织等。这其中，有的以图片形式发布，有的以二维表格形式发布，有的以数据 API 服务形式发布，比起散落在各个公告中要好多了。

数据收集的方式有线下发放问卷、从网络爬取、网络调查问卷、线下市场调查、走访、有奖征集、埋点等。在真实的数据分析中，有时候需要借助 SQL 或浏览器开发者工具，从不同的数据源获取数据，清洗整理，再将数据导入 R 环境。

2.1 从本地文件读取

利用 Base R 提供的基础函数从各类文件导入数据

2.1.1 csv 文件

小的 csv 文件，可用 Base R 提供的 read.csv() 函数读取。大型 csv 文件，可用 data.table 的 fread() 函数读取。文件 sigma.csv 存储一个 83 维的多元正态分布的协方差矩阵。

sigma <- read.csv(file = "data/sigma.csv", header = F)
sigma[1:5, 1:5]

          V1         V2         V3         V4         V5
1  0.3805689 -0.0104468  0.0000000  0.0000000  0.0000000
2 -0.0104468  0.0071204  0.0000000  0.0000000  0.0000000
3  0.0000000  0.0000000  0.3791716 -0.0016400 -0.0016488
4  0.0000000  0.0000000 -0.0016400  0.3791785 -0.0017267
5  0.0000000  0.0000000 -0.0016488 -0.0017267  0.3791763

2.1.2 xlsx 文件

readxl 读 xls 和 xlsx 文件，writexl 写 xlsx，openxlsx 读/写 xlsx 文件。

library(readxl)
cumcm2011A <- read_xls(
  path = "data/cumcm2011A附件_数据.xls",
  col_names = TRUE, sheet = "附件1", range = "A3:E322"
)
head(cumcm2011A)

# A tibble: 6 × 5
   编号 `x(m)` `y(m)` `海拔(m)` 功能区
  <dbl>  <dbl>  <dbl>     <dbl>  <dbl>
1     1     74    781         5      4
2     2   1373    731        11      4
3     3   1321   1791        28      4
4     4      0   1787         4      2
5     5   1049   2127        12      4
6     6   1647   2728         6      2

xls 文件 cumcm2011A附件_数据.xls 有多个工作表，不同的工作表有不同的数据，不同的数据占据表格的不同的位置。这可以通过函数 read_xls() 的参数 sheet 和 range 来指定工作表的名称和数据在工作表中的位置，参数 col_names 指定数据是否含有表头。

2.1.3 arrow 文件

Apache Arrow 的 R 语言接口 arrow 超出内存的大规模数据操作。比如在时空数据处理场景，数据文件往往比较大，需要在远程服务器上处理超出本地计算机内存的数据，geoarrow 包和 sfarrow 包都是应对此类需求。

2.2 从数据库中导入

从各类数据库导入数据，比如 RSQLite 等。

2.2.1 RSQLite

以 RSQLite 包内置的数据库 datasets.sqlite 为例，下面展示连接、查询和关闭连接数据库的过程。

library(DBI)
# 连接数据库
con <- dbConnect(RSQLite::SQLite(), system.file("db", "datasets.sqlite", 
        package = "RSQLite"))
# 查询数据
dbGetQuery(con, "SELECT * FROM mtcars WHERE cyl = 6 and vs = 1")

       row_names  mpg cyl  disp  hp drat    wt  qsec vs am gear carb
1 Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
2        Valiant 18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
3       Merc 280 19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
4      Merc 280C 17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4

# 参数化查询
dbGetQuery(con, 'SELECT * FROM mtcars WHERE "cyl" = :x and vs = 1',
  params = list(x = 6))

       row_names  mpg cyl  disp  hp drat    wt  qsec vs am gear carb
1 Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
2        Valiant 18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
3       Merc 280 19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
4      Merc 280C 17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4

# 断开连接
dbDisconnect(con)

2.2.2 odbc

基于开放数据库连接（Open Database Connectivity，简称 ODBC）标准，odbc 包提供了一种统一的连接各类数据库管理系统（Database Management System，简称 DBMS）的 R 接口。以连接前面的 SQLite 数据库为例，介绍 odbc 包的使用。首先安装 ODBC 驱动，接着，调用 odbc 包连接和使用数据库。在命令行中执行 odbcinst -j 可查看驱动配置文件及其位置。

在文件 odbcinst.ini 配置驱动信息，用户需根据自己的驱动安装位置调整，大致如下。

系统	驱动配置
Ubuntu	`[SQLite3] Driver = /usr/lib/x86_64-linux-gnu/odbc/libsqlite3odbc.so`
MacOS	`[SQLite3] Driver = /usr/local/lib/libsqlite3odbc.dylib`

以笔者本地的系统环境，驱动配置如下：

[SQLite3]
Driver    = /opt/homebrew/Cellar/sqliteodbc/0.99991/lib/libsqlite3odbc.dylib

在文件 odbc.ini 配置数据库连接信息（可选，可在调用函数 dbConnect() 时配置）。

[SQLite3]
Description = Test db
Driver      = SQLite3
Database    = /tmp/test
Timeout     = 2000

library(odbc)
con2 <- dbConnect(odbc::odbc(), driver = "SQLite3",
                  database = system.file("db", "datasets.sqlite", 
        package = "RSQLite") )
# 查询数据
dbGetQuery(con2, "SELECT * FROM mtcars WHERE cyl = 6 and vs = 1")

       row_names  mpg cyl  disp  hp drat    wt  qsec vs am gear carb
1 Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
2        Valiant 18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
3       Merc 280 19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
4      Merc 280C 17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4

# 参数化查询
dbGetQuery(con2, sprintf('SELECT * FROM mtcars WHERE "cyl" = %s and vs = 1', "6"))

       row_names  mpg cyl  disp  hp drat    wt  qsec vs am gear carb
1 Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
2        Valiant 18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
3       Merc 280 19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
4      Merc 280C 17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4

# 断开连接
dbDisconnect(con2)

此外，很多数据库管理系统都提供 Java 接口驱动，R 语言中也有一个 RJDBC 包，配置 Java 环境和数据库系统驱动后，可类似以上过程连接和使用数据库。

2.3 从各类网页中抓取

rvest 包从网页、网站抓取数据，再用 xml2 和 httr2 解析处理网页数据。

2.3.1 财政部官网

从中央到地方的各级政府每年都发布大量的公告，这些公告里包含大量的数据信息。比如财政部官网每年发布各类财政数据，各项收入、支出、债务以及转移支付等，下面从国务院关于中央决算的报告里，获取中央政府国债余额限额数据。

library(rvest)
# 国务院关于中央决算的报告 2012-2024 年
# https://www.mof.gov.cn/gkml/caizhengshuju/
links_vec <- c(
  `2012` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/201306/t20130629_942538.htm",
  `2013` = "https://www.mof.gov.cn/gkml/caizhengshuju/201406/t20140624_1103948.htm",
  `2014` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/201506/t20150629_1262257.htm",
  `2015` = "https://www.mof.gov.cn/gkml/caizhengshuju/201607/t20160701_2343739.htm",
  `2017` = "https://www.mof.gov.cn/gkml/caizhengshuju/201806/t20180621_2935796.htm",
  `2018` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/201906/t20190627_3286107.htm",
  `2019` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/202006/t20200622_3536392.htm",
  `2020` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/202106/t20210608_3715911.htm",
  `2021` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/202206/t20220624_3820957.htm",
  `2022` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/202307/t20230704_3894361.htm",
  `2023` = "https://www.mof.gov.cn/gkml/caizhengshuju/202407/t20240701_3938460.htm",
  `2024` = "https://www.mof.gov.cn/zhengwuxinxi/caizhengxinwen/202506/t20250630_3966883.htm"
)
# 下载网页
web_pages <- lapply(links_vec, read_html)
# 提取文本内容
extract_text <- function(x) {
  html_text(html_elements(x, "div.my_doccontent p"))
}
web_texts <- lapply(web_pages, extract_text)

rvest 包的函数 read_html 下载网页，再用函数 html_elements 选取网页中目标段落，最后，用函数 html_text 将段落转化为文本数据。有了文本数据，便可以用正则表达式提取所需数据。

# 筛选国债余额限额数据
x <- lapply(web_texts, stringr::str_extract, pattern = "(国债余额限额)(\\d+.*?)亿元", group = 2)
lapply(x, function(x) {
  x[!is.na(x)]
}) |> data.frame(check.names = F)

      2012     2013      2014      2015      2017      2018      2019      2020
1 82708.35 91208.35 100708.35 111908.35 141408.35 156908.35 175208.35 213008.35
       2021      2022      2023      2024
1 240508.35 267008.35 308608.35 352008.35

笔者曾用此方法获取历年中央决算报告中的所需数据，详见房地产相关的核心数据。

2.3.2 豆瓣排行榜

参见岁月沉淀的光影艺术，在豆瓣排行榜抓取电影Top250的数据。

2.3.3 链家二手房

2.4 从数据接口中获取

2.4.1 Github

从 Github API 接口中获取托管在 Github 上的 R 包的信息，比如点赞、关注和转发的数量。首先从 CRAN 上获得 R 包元数据信息，接着筛选出托管在 Github 上的 R 包，清理出 R 包在 Github 上的网址。

pdb <- readRDS(file = "data/cran-package-db-20231231.rds")
# 过滤出 Github 
pdb <- subset(
  x = pdb, subset = !duplicated(Package) & grepl(pattern = "github", x = BugReports),
  select = c("Package", "Maintainer", "Title", "BugReports")
)
# 掐头去尾
pdb$repo <- sub(x =  pdb$BugReports, pattern = "(http|https)://(www\\.){0,1}github\\.com/", replacement = "")
pdb$repo <- sub(x = pdb$repo, pattern = "/{1,}(issues|blob).*", replacement = "")
pdb$repo <- sub(x = pdb$repo, pattern = "/{1,}(discussions|wiki)", replacement = "")
pdb$repo <- sub(x = pdb$repo, pattern = "/$", replacement = "")

获取某代码仓库信息的 Github API 是 https://api.github.com/repos ，为了批量地访问 API ，收集想要的数据，将数据请求、结果整理的过程打包成一个函数。

github_stats <- function(repo) {
  url <- paste("https://api.github.com/repos", repo, sep = "/")
  # 最多允许失败 5 次，每失败一次休息 5s
  req <- xfun::retry(curl::curl_fetch_memory, url = url, .times = 5, .pause = 5)
  x <- jsonlite::fromJSON(rawToChar(req$content))
  # 爬失败的标记一下
  if(is.null(x$stargazers_count)) x$stargazers_count <- x$subscribers_count <- x$forks_count <- -1
  # 爬一个休息 1s
  Sys.sleep(1)
  data.frame(
    repo = repo,
    # 点赞 仓库上 star 的人数
    stargazers_count = x$stargazers_count,
    # 关注 仓库上 watch 的人数
    subscribers_count = x$subscribers_count,
    # 转发 仓库上 fork 的人数
    forks_count = x$forks_count
  )
}

下面测试一下这段代码，获取代码仓库 yihui/knitr 的点赞、关注和转发的人数。

# 测试代码
github_stats(repo = "yihui/knitr")

         repo stargazers_count subscribers_count forks_count
1 yihui/knitr               -1                -1          -1

理论上，使用函数 lapply() 遍历所有 R 包可得所需数据，将数据收集函数应用到每一个 R 包上再合并结果，即如下操作。

# 合并数据
gh_repo_db <- data.table::rbindlist(lapply(pdb$repo, github_stats))

实际上，在没有访问令牌的情况下，Github API 的访问次数是有限制的，只有 60 次（一段时间内）。首先在 Github 开发者设置中申请一个应用，获得应用名称（appname）、客户端 ID（key）和密钥（secret），下面借助 httr 包配置 OAuth 凭证。

library(httr)
# Github API Oauth2
oauth_endpoints("github")
# 应用名称（appname）、客户端 ID（key）和密钥（secret）
myapp <- oauth_app(
  appname = "Application Name", key = "Client ID",
  secret = "Client Secrets"
)
# 获取 OAuth 凭证
github_token <- oauth2.0_token(oauth_endpoints("github"), myapp)
# 使用 API
gtoken <- config(token = github_token)

修改函数 github_stats() 中请求 Github API 的一行代码，发送带密钥的 GET 请求。

req <- xfun::retry(GET, url = url, config = gtoken, .times = 5, .pause = 5)

此外，请求难免出现意外，按照上面的方式，一旦报错，数据都将丢失。因此，要预先准备存储空间，每获取一条数据就存进去，如果报错了，就打个标记。

# 准备存储数据
gh_repo_db <- data.frame(
  repo = pdb$repo, stargazers_count = rep(-1, length(pdb$repo)),
  subscribers_count = rep(-1, length(pdb$repo)),
  forks_count = rep(-1, length(pdb$repo))
)
# 不断更新数据
while (any(gh_repo_db$stargazers_count == -1)) {
  tmp <- gh_repo_db[gh_repo_db$stargazers_count == -1, ]
  for (repo in tmp$repo) {
    gh_repo_db[gh_repo_db$repo == repo, ] <- github_stats(repo = repo)
  }
  if(repo == tmp$repo[length(tmp$repo)]) break
}

最后，将收集整理好的数据保存到磁盘上，下面按点赞数量给 R 包排序，篇幅所限，仅展示前 20。

gh_repo_db <- readRDS(file = "data/gh-repo-db-2023.rds")
gh_repo_db <- gh_repo_db[!duplicated(gh_repo_db$repo),]
gh_repo_db <- gh_repo_db[order(gh_repo_db$stargazers_count, decreasing = T),] 
head(gh_repo_db, 20)

                      repo stargazers_count subscribers_count forks_count
8434          dmlc/xgboost            25266               909        8707
5553      facebook/prophet            17415               425        4474
4307         mlflow/mlflow            16365               292        3793
3807    Microsoft/LightGBM            15821               437        3798
265           apache/arrow            13080               356        3220
3080           h2oai/h2o-3             6624               384        2016
2790     tidyverse/ggplot2             6210               308        2028
3430 interpretml/interpret             5894               141         706
6921         rstudio/shiny             5180               339        1818
4317         mlpack/mlpack             4668               185        1577
1754       tidyverse/dplyr             4612               246        2131
640       rstudio/bookdown             3565               122        1263
1430 Rdatatable/data.table             3437               170         977
6316     rstudio/rmarkdown             2758               146         977
2269          wesm/feather             2708                97         174
5324       plotly/plotly.R             2467               117         628
5084   thomasp85/patchwork             2344                49         159
1389        r-lib/devtools             2340               120         760
3658           yihui/knitr             2326               115         877
6868      satijalab/seurat             2034                75         867

将发布在 Github 上的受欢迎的 R 包列出来了，方便读者选用，也看到一些有意思的结果。

机器学习相关的 R 包靠在最前面，实际上，它们（占十之七八）多是对应软件的 R 语言接口，点赞的数目应当算上其它语言接口的贡献。
在机器学习之后，依次是数据可视化（ggplot2、shiny、plotly.R、patchwork）、数据操作（dplyr、data.table、feather）和可重复性计算（bookdown、rmarkdown、knitr）、R 包开发（devtools）和生物信息（seurat）。

最后，简要说明数据的情况：以上观察结果是基于 CRAN 在 2023-12-31 发布的 R 包元数据，8475 个 R 包在 Github 托管源代码，这些 R 包的点赞、关注和转发数据是在 2024-01-30 爬取的。其中，共有 29 个 R 包不按规矩填写、改名字、换地方、甚至删库了，这些 R 包是可忽略的。当然，也存在一些 R 包并未托管在 Github 上，但质量不错，比如 glmnet 包、colorspace 包、fGarch 包、bnlearn 包等，应当是少量的。

2.4.2 中国地震台网

中国地震台网可以想象后台有一个数据库，在页面的小窗口中输入查询条件，转化为某种 SQL 语句，传递给数据库管理系统，执行查询语句，返回查询结果，即数据。类似地，美国地质调查局提供一些选项窗口，可供选择数据范围，直接下载 CSV 或 XLS 文件。笔者曾下载了2012-2022 年的地震数据，并用各类可视化手段探索数据，发现地下很深的位置往往发生大地震，而大地震不一定发生在地下很深的位置。详见地震越来越频繁了吗？

2.4.3 美国人口调查局

美国人口调查局每年会发布大量国家社会经济发展水平相关的数据，通过访问令牌可以批量下载数据。用户需要先注册账号，获取使用 API 接口的访问令牌，可以想象后台不仅有一个数据库，在此之上，还有一层数据鉴权。tidycensus 包将数据下载并整理成 tidyverse 风格数据集。

library(tidycensus)
library(sf)
options(tigris_use_cache = TRUE)
# 下载北卡州郡级人口调查数据
nc_income_race_county <- get_acs(
  state = "NC", # 北卡州
  geography = "county", # 郡级
  # 三个数据指标：家庭年收入，白人数量，总人数
  variables = c("B19013_001", "B02001_001", "B02001_002"),
  key = Sys.getenv("CENSUS_API_KEY"), # 设置访问令牌
  geometry = TRUE, # 调 tigris 包下载北卡州郡级地图数据 
  progress_bar = F, # 不显示下载进度
  output = "wide", # 输出宽格式数据
  year = 2023,     # 2019-2023 年度
  survey = "acs5", # 5 年
  moe_level = 90   # 置信水平 90%
)
nc_income_race_county

Simple feature collection with 100 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -84.32187 ymin: 33.84232 xmax: -75.46062 ymax: 36.58812
Geodetic CRS:  NAD83
First 10 features:
   GEOID                               NAME B19013_001E B19013_001M B02001_001E
1  37133      Onslow County, North Carolina       64568        1873      208537
2  37009        Ashe County, North Carolina       50827        3707       26831
3  37169      Stokes County, North Carolina       60039        2412       44889
4  37053   Currituck County, North Carolina       91548        8181       29612
5  37173       Swain County, North Carolina       55429        6686       14065
6  37131 Northampton County, North Carolina       47935        3311       17212
7  37153    Richmond County, North Carolina       43626        5421       42818
8  37051  Cumberland County, North Carolina       58780         855      336749
9  37085     Harnett County, North Carolina       69012        3266      136503
10 37199      Yancey County, North Carolina       54961        5646       18676
   B02001_001M B02001_002E B02001_002M                       geometry
1           NA      144394        1126 MULTIPOLYGON (((-77.17131 3...
2           NA       25215         344 MULTIPOLYGON (((-81.74065 3...
3           NA       40992         304 MULTIPOLYGON (((-80.4502 36...
4           NA       25360         338 MULTIPOLYGON (((-76.3133 36...
5           NA        8601         125 MULTIPOLYGON (((-83.94939 3...
6           NA        6827          90 MULTIPOLYGON (((-77.89977 3...
7           NA       24360         363 MULTIPOLYGON (((-80.07567 3...
8           NA      147878        1159 MULTIPOLYGON (((-79.11285 3...
9           NA       85935         904 MULTIPOLYGON (((-79.22186 3...
10          NA       17268         127 MULTIPOLYGON (((-82.50538 3...

下载了2019-2023年度美国北卡罗来纳州各个郡的人口调查统计数据，分别是家庭年收入、白人数量和总人数，同时下载北卡州郡级行政区划地图数据。笔者曾用此方法获取数据，可视化北卡罗来纳州家年年收入和白人占比的关系，详见地区分布图及其应用。

2.4.4 世界银行

世界银行和国际货币基金组织都公开一些经济数据，wbstats 包和 WDI 包封装世界银行提供的数据接口 REST API ，下面通过 wbstats 包下载2024 年全球各国的预期寿命、人均 GDP和人口总数数据。

library(wbstats)
# 从世界银行数据库获取宏观经济数据
wb_world_indicator <- wbstats::wb_data(
  indicator = c(
    "SP.DYN.LE00.IN", # 预期寿命
    "NY.GDP.PCAP.CD", # 人均 GDP
    "SP.POP.TOTL" # 人口总数
  ),
  country = "countries_only",
  return_wide = FALSE,
  start_date = 2010, end_date = 2023
)
wb_world_indicator

# A tibble: 9,114 × 10
   indicator_id   indicator     iso2c iso3c country  date value unit  obs_status
   <chr>          <chr>         <chr> <chr> <chr>   <dbl> <dbl> <chr> <chr>     
 1 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2023  66.0 <NA>  <NA>      
 2 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2022  65.6 <NA>  <NA>      
 3 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2021  60.4 <NA>  <NA>      
 4 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2020  61.5 <NA>  <NA>      
 5 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2019  62.9 <NA>  <NA>      
 6 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2018  62.4 <NA>  <NA>      
 7 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2017  62.4 <NA>  <NA>      
 8 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2016  62.6 <NA>  <NA>      
 9 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2015  62.3 <NA>  <NA>      
10 SP.DYN.LE00.IN Life expecta… AF    AFG   Afghan…  2014  62.3 <NA>  <NA>      
# ℹ 9,104 more rows
# ℹ 1 more variable: last_updated <date>

笔者曾经从这个数据源去复现汉斯·罗琳的动画，详见Gapminder：关注差异、理解变化。