• Posts tagged "数据"

Blog Archives

用R语言获取全球新冠疫情数据

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统  计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-coivd19-data-nCov2019

前言

新冠长期以来影响着我们的工作和生活,近期在北京和上海两个超大城市,还在不停的变种和传播,居家办公,每天核酸,都已经变成了生活的常态。

对新冠数据进行数据分析,就需要获取新冠的数据。中国的团队开发的R包nCov2019,可以方便的获得新冠全量数据,帮助我们提高数据获取的效率。

目录

  1. 开源新冠疫情数据包 nCov2019
  2. 安装nCov2019包
  3. 查看数据集
  4. 疫情可视化:世界地图
  5. 疫情可视化:多种分析图表
  6. 疫情可视化:shiny控制台

1. 获取新冠疫情数据

当前有很多的新冠疫情开源的数据源,为了简化数据获取难度、并分析方便,我找到一个R语言新冠数据开源项目 YuLab-SMU/nCov2019。nCov2019包,收集了从全球新冠疫情的数据,从2020年01月22日开始到当前的每日新冠的数据,支持R语言的API的,我们直接使用这个包就可以获得全球数据,并在github上的开源。

nCov2019项目,为了方便获取有关冠状病毒爆发的流行病学数据,除了详细的基础统计数据外,它还包括有关疫苗开发和候选疗法的信息。 并设计了用于地理地图可视化的函数 plot() 并提供了一个交互式shiny应用程序。 这些功能有助于告知公众和研究这种病毒和类似病毒如何在人口稠密的国家传播。

nCov2019项目,是由中国的南方医科大学基础医学院生物信息学系研发,研发团队:https://yulab-smu.top/members/。给中国的团伙点 “赞”。

2.安装nCov2019包

首先,安装nCov2019包,如果一切顺利,一条命令就完成了。


# 如果没有remotes包,请先安装remotes包,这样才能从github上下载代码。
# > install.packages("remotes")

# 安装YuLab-SMU/nCov2019包
> remotes::install_github("YuLab-SMU/nCov2019")
Using github PAT from envvar GITHUB_PAT
Downloading GitHub repo YuLab-SMU/nCov2019@HEAD

小插曲:如果在安装时,出现 GITHUB_PAT错误时。

Using github PAT from envvar GITHUB_PAT
Error: Failed to install 'unknown package' from GitHub:
  HTTP error 401.
  Bad credentials

  Rate limit remaining: 27/60
  Rate limit reset at: 2021-09-30 19:17:49 UTC

说明你的GITHUB没有设置 开通Personal access tokens的权限,请访问打开https://github.com/settings/tokens 点击 generate new taken 按钮,创建一个新的token,复制这个新生成的token。之后,切回到R语言开发环境,使用R语言命令,会打开一个新文件,并输入新生新的token。

> usethis::edit_r_environ()
* Modify 'C:/Users/bsspi/Documents/.Renviron'
* Restart R for changes to take effect

文件内容:请替换自己的token。

GITHUB_PAT=ghp_a4e7l0KjhVVzX5J86PjVcFS7A7u8st1rvQIN

设置完GITHUB_PAT后,保存再重启R语言环境,再执行安装命令就可以正常下载和安装了。

> remotes::install_github("YuLab-SMU/nCov2019")

加载nCov2019包,再加载数据

> library(nCov2019)
> x <- query()
Querying the latest data...
last update: 2022-07-05 
Querying the global data...
Gloabl total  555324897  cases; and  6362663  deaths
Gloabl total affect country or areas: 230
Gloabl total recovered cases: 165929
last update: 2022-07-05 
Querying the historical data...
Querying the vaccine data...
Total Candidates Programs : 51 
Querying the therapeutics data...
Total Candidates Programs : 84 
Query finish, each time you can launch query() to reflash the data

#查看x对象
> x
$latest
last update: 2022-07-05 

$global
$updated
[1] "2022-07-05"

$cases
[1] 555324897

$todayCases
[1] 137897

$deaths
[1] 6362663

$todayDeaths
[1] 261

$recovered
[1] 530092910

$todayRecovered
[1] 165929

$active
[1] 18869324

$critical
[1] 38501

$casesPerOneMillion
[1] 71243

$deathsPerOneMillion
[1] 816.3

$tests
[1] 6480187084

$testsPerOneMillion
[1] 816521

$population
[1] 7936338827

$oneCasePerPeople
[1] 0

$oneDeathPerPeople
[1] 0

$oneTestPerPeople
[1] 0

$activePerOneMillion
[1] 2377.59

$recoveredPerOneMillion
[1] 66793.13

$criticalPerOneMillion
[1] 4.85

$affectedCountries
[1] 230

attr(,"class")
[1] "global_summary"
attr(,"row.names")
[1] 1

$historical
last update: 2022-07-04 

$vaccine
Total Candidates Programs : 51 

$therapeutics
Total Candidates Programs : 84 

由下载数据比较慢,所以我们可以下载完成后,先在本地存一份

# 设置工作路径
> setwd("C:/work/R/covid19")

# 本地存储
# > saveRDS(x,"x.rds")

# 本地读取
# > x<-readRDS("x.rds")

3. 查看数据集

查看 x 对象变量,x对象是一个list类型,包括了全部的新冠数据集,按使用方便主要分另5个数据集进行存储。我们可以用names()函数,找到这5个数据集。

> names(x)
[1] "latest"       "global"       "historical"   "vaccine"      "therapeutics"
  • global,全球总体汇总统计
  • latest,全球所有国家的最新一天数据
  • historical,全球所有国家的历史数据
  • vaccine,目前疫苗研发进展
  • therapeutics,目前的治疗学发展进展

3.1 全球总体汇总统计 x$global
查看 x$global 的汇总信息,今天日期为20220705,为最后更新日期。累计总病例为 555324897 例 其中死亡 6362663 例,设计全球国家230个,治愈165929例。

> global<-x$global
> summary(global)
Gloabl total  555324897  cases; and  6362663  deaths
Gloabl total affect country or areas: 230
Gloabl total recovered cases: 165929
last update: 2022-07-05

# 全球致死率
> global$deaths/global$cases
[1] 0.01145755

3.2 全球所有国家的最新一天数据 x$latest
查看 x$latest 的截止到最新的累计数据,打印10条数据。更新数据来自:https://www.worldometers.info/coronavirus/

> last<-x$latest
> head(last["Global"],10)
     country    cases  deaths recovered  active todayCases todayDeaths todayRecovered population      tests    updated
1    Germany 28542484  141397  26886400 1514687     147489         102         143600   84320113  122332384 2022-07-05
2      Italy 18805756  168604  17617973 1019179      36282          59          27590   60283926  227930923 2022-07-05
3     Brazil 32536302  672101  30967114  897087      33833          84          60539  215583725   63776166 2022-07-05
4  Australia  8254785   10039   7979050  265696      29661          25          17866   26092113   74123177 2022-07-05
5     France 31452173  149726  29774233 1528214      24418          97          73899   65562570  271490188 2022-07-05
6     Taiwan  3893612    7025   2996511  890076      23087          69          66067   23903044   21241708 2022-07-05
7      Japan  9398126   31316   9179820  186990      22154           8          13900  125701774   56860191 2022-07-05
8     Israel  4391284   10984   4264371  115929      19378           7              0    9326000   41373364 2022-07-05
9      India 43532788  525242  42891933  115613      14224          19          12456 1407162089  863999907 2022-07-05
10       USA 89567321 1043372  85286101 3237848      13317          38         113249  334887583 1052987802 2022-07-05

选择美国、印度、中国的汇总数据。

> last[c("USA","India","China")]
   country    cases  deaths recovered  active todayCases todayDeaths todayRecovered population      tests    updated
9    India 43532788  525242  42891933  115613      14224          19          12456 1407162089  863999907 2022-07-05
10     USA 89567321 1043372  85286101 3237848      13317          38         113249  334887583 1052987802 2022-07-05
76   China   225923    5226    220165     532         72           0             50 1439323776  160000000 2022-07-05

查看中国在2022-07-05日的明细数据。

> last$detail[which(last$detail$country=="China"),]
      updated country countryInfo._id countryInfo.iso2 countryInfo.iso3 countryInfo.lat countryInfo.long
76 2022-07-05   China             156               CN              CHN              35              105
                             countryInfo.flag  cases todayCases deaths todayDeaths recovered todayRecovered active critical
76 https://disease.sh/assets/img/flags/cn.png 225923         72   5226           0    220165             50    532        0
   casesPerOneMillion deathsPerOneMillion     tests testsPerOneMillion population continent oneCasePerPeople oneDeathPerPeople
76                157                   4 160000000             111163 1439323776      Asia             6371            275416
   oneTestPerPeople activePerOneMillion recoveredPerOneMillion criticalPerOneMillion
76                9                0.37                 152.96                     0

3.3 全球所有国家的历史数据x$historical
查看历史数据,数据源来自:https://coronavirus.jhu.edu/map.html

# 创建变量
> hist<-x$historical

# 查看中国的历史数据前10行
> head(hist["China"],10)
     country       date cases deaths recovered
39     China 2020-01-22   548     17        28
238    China 2020-01-23   643     18        30
437    China 2020-01-24   920     26        36
636    China 2020-01-25  1406     42        39
835    China 2020-01-26  2075     56        49
1034   China 2020-01-27  2877     82        58
1233   China 2020-01-28  5509    131       101
1432   China 2020-01-29  6087    133       120
1631   China 2020-01-30  8141    171       135
1830   China 2020-01-31  9802    213       214

# 查看中国的历史数据后10行
> tail(hist["China"],10)
       country       date   cases deaths recovered
176154   China 2022-06-25 2124209  14624         0
176353   China 2022-06-26 2126233  14624         0
176552   China 2022-06-27 2128225  14624         0
176751   China 2022-06-28 2130033  14625         0
176950   China 2022-06-29 2132190  14625         0
177149   China 2022-06-30 2134718  14627         0
177348   China 2022-07-01 2137169  14628         0
177547   China 2022-07-02 2139919  14631         0
177746   China 2022-07-03 2142216  14633         0
177945   China 2022-07-04 2144566  14633         0

查看中国北京的数据

> head(hist['China','beijing'])
    country province       date cases deaths recovered
6     China  beijing 2020-01-22    14      0         0
95    China  beijing 2020-01-23    22      0         0
184   China  beijing 2020-01-24    36      0         1
273   China  beijing 2020-01-25    41      0         2
362   China  beijing 2020-01-26    68      0         2
451   China  beijing 2020-01-27    80      1         2

可视化北京疫情数据情况


> library(ggplot2)
> library(reshape2)
> beijing_df<-hist['China','beijing']
> beijing<-melt(beijing_df[,c("date", "cases" ,"deaths", "recovered")],id.vars = c("date"))
> 
> g<-ggplot(data = beijing, mapping = aes(x=date,y=value,colour=variable))
> g<-g+geom_line() + geom_point()                                   #绘制线图和点图
> g<-g+scale_shape_manual(values = c(21,23))                        #自定义点形状
> g<-g+scale_y_log10()
> g<-g+ggtitle("北京疫情统计")+xlab("日期")+ylab("Log(人数)")
> g

上图中,红色线为累计感染新冠人数,绿色线为累计致死人数,蓝色线为治愈人数,从2021-08-05 数据就停了,因为出现了蓝色线直接到0的情况。

3.4 目前疫苗研发进展 x$vaccine

nCov2019数据集,还收集了疫苗研发进展等相关数据,数据源来自:https://www.raps.org/news-and-articles/news-articles/2020/3/covid-19-vaccine-tracker

用户可以检查疫苗的开发状态,使用summary() 来查看疫苗的试验阶段,和候选疫苗数量。

# 查看疫苗数据
> vac <- x$vaccine

# 尝试阶段,和候选疫苗
> summary(vac)
         phase candidates
1      Phase 3         10
2    Phase 2/3          3
3      Phase 2          2
4    Phase 1/2          9
5      Phase 1         13
6 Pre-clinical         14

查看每种疫苗的摘要信息,例如机制、试验阶段、 机构等。

> head(vac["all"])
   id candidate                                                                mechanism                 sponsors trialPhase
1 id1    BNT162                                                       mRNA-based vaccine         Pfizer, BioNTech    Phase 3
2 id2 mRNA-1273                                                       mRNA-based vaccine                  Moderna    Phase 3
3 id3  Ad5-nCoV                           Recombinant vaccine (adenovirus type 5 vector)        CanSino Biologics    Phase 3
4 id4   AZD1222 Replication-deficient viral vector vaccine (adenovirus from chimpanzees) The University of Oxford    Phase 3
5 id5 CoronaVac                        Inactivated vaccine (formalin with alum adjuvant)                  Sinovac    Phase 3
6 id6   Covaxin                                                      Inactivated vaccine           Bharat Biotech    Phase 3
                                             institutions
1 Multiple study sites in Europe, North America and China
2  Kaiser Permanente Washington Health Research Institute
3                                         Tongji Hospital
4                          The University of Oxford, 
5              Sinovac Research and Development Co., Ltd.
6                     

在中国,我们主要使用的 科兴生物公司的疫苗。科兴控股生物技术有限公司(简称科兴生物,英语:Sinovac Biotech)是中国的生物技术和疫苗生产公司,总部设于北京市海淀区。科兴生物成立于1999年,曾在美国上市,公司通过全资子公司科兴控股(香港)有限公司,拥有北京科兴生物制品有限公司、科兴(大连)疫苗技术有限公司、北京科兴中维生物技术有限公司和北京科兴中益生物医药有限公司四家企业。

科兴2019冠状病毒病疫苗(英语:Sinovac COVID-19 Vaccine,俗称科兴疫苗,商品名克尔来福(英语:CoronaVac)是中国大陆生物制药公司科兴生物开发的2019冠状病毒病疫苗,属于灭活疫苗。2021年6月1日,世卫组织宣布将克尔来福列入“紧急使用清单”。


# 疫苗名称 CoronaVac
> vac["all"][which(vac["all"]$candidate=='CoronaVac'),]
   id candidate                                         mechanism sponsors trialPhase                               institutions
5 id5 CoronaVac Inactivated vaccine (formalin with alum adjuvant)  Sinovac    Phase 3 Sinovac Research and Development Co., Ltd.

查看id5的科兴疫苗的详细背景信息。


> vac[ID="id5"]
[1] "Background: CoronaVac (formerly PiCoVacc) is a formalin-inactivated 
and alum-adjuvanted candidate vaccine. Results from animal studies showed 
“partial or complete protection in macaques” exposed to SARS-CoV-2, 
according to a paper published  in Science.  Study Design: A 
Phase 1/2 trial of 743 healthy volunteers (18-59 years old) who received 
two different dosages of the vaccine or placebo is active but not 
recruiting. A Phase 1 trial of 143 participants (NCT04352608) and a Phase 2 
trial of 600 participants (NCT04383574) are both active but not recruiting. 
Sinovac said a Phase 3 trial in collaboration with Instituto Butantan in 
Brazil is underway (NCT04456595), and the company plans to enroll around 
9,000 patients in the healthcare industry. Trials also are underway in 
Turkey (NCT04582344) and in Indonesia (NCT04508075).  Outcomes: Results 
from the Phase 1/2 trials published in The Lancet Infectious Diseases 
indicate the vaccine has good safety and immunogenicity, with 
seroconversion occurring in 92.4% of participants receiving the 3 μg dose 
on a 0-14 day schedule and 97.4% of individuals receiving the same dose on 
a 0-28 day schedule.  Status:  Representatives from Sinovac told Reuters 
that the vaccine appeared to be safe in older trial participants, and did 
not cause any severe side effects. Preliminary results from the Instituto 
Butantan trial announced by the company indicate CoronaVac is safe so far, 
with no serious adverse events reported. The trial in Brazil was briefly 
suspended due to a patient death, but the trial has since resumed."

Google翻译后:“背景:CoronaVac(前身为 PiCoVacc)是一种福尔马林灭活和明矾佐剂的候选疫苗。根据一篇论文,动物研究结果显示暴露于 SARS-CoV-2 的猕猴“部分或完全保护”发表在《科学》杂志上。研究设计:接受两种不同剂量疫苗或安慰剂的 743 名健康志愿者(18-59 岁)的 1/2 期试验处于活跃状态,但未招募。143 名参与者的 1 期试验 (NCT04352608) 600 名参与者的第二阶段试验 (NCT04383574) 均处于活跃状态,但未招募。科兴生物表示,与巴西 Instituto Butantan 合作的第三阶段试验正在进行中 (NCT04456595),该公司计划在医疗保健行业招募约 9,000 名患者. 土耳其 (NCT04582344) 和印度尼西亚 (NCT04508075) 的试验也在进行中。 结果: 发表在《柳叶刀传染病》上的 1/2 期试验结果表明,该疫苗具有良好的安全性和免疫原性,具有良好的安全性和免疫原性。 92.4% 的参与者在 0-14 天接受 3 μg 剂量,97.4% 的个体在 0-28 天接受相同剂量。状态:科兴生物的代表告诉路透社,该疫苗在年长的试验参与者中似乎是安全的,并且没有引起任何严重的副作用。该公司宣布的 Instituto Butantan 试验的初步结果表明 CoronaVac 到目前为止是安全的,没有报告严重的不良事件。由于患者死亡,巴西的试验曾短暂暂停,但此后又恢复了试验。”

3.5 目前的治疗学发展进展x$therapeutics

nCov2019包,治疗学发展进展的相关数据,数据源来自:https://aps.org/news-and-articles/news-articles/2020/3/covid-19-therapeutics-tracker

用户可以查看药品的开发状态,使用summary() 来查看治疗的药品试验阶段和候选药品数量。

#查看治疗数据
> thera<- x$therapeutics
> summary(thera)
                                    phase candidates
1                                 Phase 3         13
2                             Phase 2/3/4          3
3                               Phase 2/3         28
4                               Phase 1/2          1
5                                 Phase 2         15
6                               Phase 3/4          2
7    No longer being studied for COVID-19          4
8                                 Various          1
9                                 Phase 1          4
10                             Phase 2b/3          2
11 No longer being developed for COVID-19          1
12                            Phase 1/2/3          1
13                              Phase 1/4          1
14                            Phase 1b/2a          1
15                                Phase 4          1
16                              Phase 2/2          1
17                               Phase 1b          4
18                              Phase 2/4          1

查看所有药品

> head(thera["All"])
   id       medicationClass                                                        tradeName       developerResearcher
1 id1             Antiviral Molnupiravir (Lagevrio, formerly known as MK-4482 and EIDD-2801) Ridgeback Biotherapeutics
2 id2   Monoclonal antibody                             Evusheld (tixagevimab and cilgavimab               AstraZeneca
3 id3   Monoclonal antibody                                  Regkirona (regdanvimab, CT-P59)                 Celltrion
4 id4 IL-6 receptor agonist                                  Actemra/RoActemra (tocilizumab)                     Roche
5 id5   Monoclonal antibody    Amubarvimab and romlusevimab (formerly BRII-196 and BRII-198)  Brii Biosciences Limited
6 id6         Anticoagulant                                             Heparin (UF and LMW)                     NHLBI
                   sponsors  trialPhase lastUpdate
1 Ridgeback Biotherapeutics     Phase 3 2020-12-10
2               AstraZeneca     Phase 3 2020-12-10
3                 Celltrion     Phase 3 2020-12-10
4                   Various     Phase 3 2020-12-10
5                     NIAID     Phase 3 2020-12-10
6      Operation Warp Speed Phase 2/3/4 2020-12-10

查看辉瑞公司开发了一种名为 Paxlovid 的 COVID-19 抗病毒治疗组合。 Paxlovid 是 nirmatrelvir(原 PF-07321332)(一种蛋白酶抑制剂,可阻断 SARS-CoV-2 中的一种酶并阻止病毒复制)和 HIV 抗病毒药物利托那韦的组合。


# 药品名称 Oral antiviral
> thera["All"][which(thera["All"]$medicationClass=='Oral antiviral'),]
     id medicationClass                tradeName                developerResearcher                           sponsors trialPhase
30 id30  Oral antiviral TEMPOL (4-Hydroxy-TEMPO) Adamis Pharmaceuticals Corporation Adamis Pharmaceuticals Corporation  Phase 2/3
   lastUpdate
30 2020-10-08

查看对本药品的背景介绍。


> thera[ID="id30"] 
[1] "Background: TEMPOL is an oral antiviral drug being evaluated as a 
potential at-home therapeutic for COVID-19. Researchers believe TEMPOL can 
help to break up clusters of iron and sulfur cells that SARS-CoV-2 uses to 
replicate. Additionally, researchers noted TEMPOL has anti-cytokine and 
anti-inflammatory properties.Trial: A Phase 2/3 trial is currently 
recruiting up to 248 participants to test whether TEMPOL can limit 
hospitalization in patients with an early SARS-CoV-2 infection 
(NCT04729595). At least 50 individuals in the study enrolled will have 
comorbidities, according to the ClinicalTrials.gov listing."

Google翻译:TEMPOL 是一种口服抗病毒药物,被评估为 COVID-19 的潜在家庭治疗药物。研究人员认为,TEMPOL 可以帮助分解 SARS-CoV-2 用来复制的铁和硫细胞簇。 此外,研究人员指出 TEMPOL 具有抗细胞因子和抗炎特性。试验:2/3 期试验目前正在招募多达 248 名参与者,以测试 TEMPOL 是否可以限制早期 SARS-CoV-2 感染患者的住院治疗 (NCT04729595 ). 根据 ClinicalTrials.gov 的列表,参与研究的至少 50 人将患有合并症。”

4. 疫情数据可视化展示

nCov2109包,不仅可以让我们免费的获取新冠相关的数据,还支持一些可视化的展示,方便我们理解和使用数据。

全球累计确诊病例,按国家进行可视化。

> plot(last)

全球累计新冠检测数据,按国家进行可视化。

> plot(last, type="tests",palette="Green")

比较不同国家每天新增确诊病例的数量。

# 加载画图包
> library(ggplot2)
> library(dplyr)

# 构建数据
> tmp <- hist["global"] %>%
+   group_by(country) %>%
+   arrange(country,date) %>%
+   mutate(diff = cases - lag(cases, default =  first(cases))) %>%
+   filter(country %in% c("Australia", "Japan", "Italy", "Germany",  "China")) 

# 画图
> ggplot(tmp,aes(date, log(diff+1), color=country)) + geom_line() +
+   labs(y="Log2(daily increase cases)") + 
+   theme(axis.text = element_text(angle = 15, hjust = 1)) +
+   scale_x_date(date_labels = "%Y-%m-%d") + 
+   theme_minimal()

指定日期,用历史数据绘制过去时间的爆发图。

> plot(hist, region="Global" ,date = "2020-08-01", type="cases")

动画效果图:画从出2022-03-01 到 2022-05-01 的数据,输出为 GIF 动画。

> from = "2020-03-01"
> to = "2020-04-01"
> plot(hist, from = from, to=to)

5. 疫情可视化:多种分析图表

nCov2019包,同时提供了其他维度的可视化分析功能。

中国累计确诊可视化。


> china <- hist['China']
> china <- china[order(china$cases), ]
> ggplot(china, 
+        aes(date, cases)) +
+   geom_col(fill = 'firebrick') + 
+   theme_minimal(base_size = 14) +
+   xlab(NULL) + ylab(NULL) + 
+   scale_x_date(date_labels = "%Y/%m/%d") +
+   labs(caption = paste("accessed date:", max(china$date)))

画图展示8个国家的新冠确珍人数每日增速,并进行线性拟合。

> library(dplyr)
> library(magrittr)
> library(ggrepel)

# 指定国家 
> country_list =  c("Italy","Brazil","Japan","China","USA","Mexico","Russia","India","Thailand")

# 数据处理 
> hist[country_list]  %>%
+     subset( date > as.Date("2020-10-01") ) %>%
+     group_by(country) %>%
+     arrange(country,date) %>%
+     mutate(increase = cases - lag(cases, default =  first(cases))) -> df

# 画图 
> ggplot(df, aes(x=date, y=increase, color=country  ))+
+     geom_smooth() + 
+     geom_label_repel(aes(label = paste(country,increase)), 
+                      data = df[df$date == max(df$date), ], hjust = 1) + 
+     labs(x=NULL,y=NULL)+ 
+     theme_bw() + theme(legend.position = 'none') 

从图中可以直观看到,截止到2022年7月,当前选中的8个国家,阳性病例每日确诊都已经过了高速增长阶段,增速已明显下降。

画出印度新冠确诊和治愈的趋势。


> library('tidyr')
> country<-"India"
> india<-hist[country]
> india <- gather(india, curve, count, -date, -country)
> 
> ggplot(india, aes(date, count, color = curve)) + geom_point() + geom_line() + 
+   labs(x=NULL,y=NULL,title=paste("Trend of cases, recovered and deaths in", country)) +
+   scale_color_manual(values=c("#f39c12", "#dd4b39", "#00a65a")) +
+   theme_bw() +   
+   geom_label_repel(aes(label = paste(curve,count)), 
+                    data = india[india$date == max(india$date), ], hjust = 1) + 
+   theme(legend.position = "none", axis.text = element_text(angle = 15, hjust = 1)) +
+   scale_x_date(date_labels = "%Y-%m-%d")

通过热力图画出,199个国家,新冠确诊病例影响程度。

> all <- hist["global"]
> all <- all[all$cases > 0,]
> length(unique(all$country))
[1] 199

> all <- subset(all,date <= as.Date("2020-3-19"))
> max_time <- max(all$date)
> min_time <- max_time - 7
> all <-  all[all$date >= min_time,]
> all2 <- all[all$date == max(all$date,na.rm = TRUE),]

> all$country <- factor(all$country, levels=unique(all2$country[order(all2$cases)]))
> breaks = c(0,1000, 10000, 100000, 10000000)

> ggplot(all, aes(date, country)) +
+   geom_tile(aes(fill = cases), color = 'black') +
+   scale_fill_viridis_c(trans = 'log', breaks = breaks, labels = breaks) +
+   xlab(NULL) + ylab(NULL) +
+   scale_x_date(date_labels = "%Y-%m-%d") + theme_minimal()

从图中对比来看,中国是防疫最高的国家。

6. 疫情可视化:shiny控制台

nCov2019包,还提供了shiny的dashboard功能,面向最终使用者,生成动态可交互的效果,可以直接用于功能演示,体现R语言的强大的原型功能。

# 打开shiny 控制台
> dashboard()

本文详细地介绍了nCov2019包的使用,我们可以通过nCov2019包来方便地获取全球的新冠的数据,并且可以按国家按地区按时间,进行不同数据维度的提取。同时nCov2019包,还提供了可视化分析和可视化的shiny控制台,能帮助用户直观的理解数据,而不需要写大量程序。

文本完整的代码,已经上传github,可以自由下载使用:https://github.com/bsspirit/infect/blob/main/code/getdata.r

转载请注明出处:
http://blog.fens.me/r-coivd19-data-nCov2019

打赏作者

2020 Microsoft Ignite The Tour ShenZhen

跨界知识聚会系列文章,“知识是用来分享和传承的”,各种会议、论坛、沙龙都是分享知识的绝佳场所。我也有幸作为演讲嘉宾参加了一些国内的大型会议,向大家展示我所做的一些成果。从听众到演讲感觉是不一样的,把知识分享出来,你才能收获更多。

关于作者

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/meeting-mvp-20200114

前言

2020年,微软在全球30个城市开展Microsoft Ignite The Tour活动,Microsoft专家和技术社区将为您带来最好的Microsoft Ignite活动体验,让您了解构建解决方案以及迁移和管理基础结构的新方法,并与当地行业领导者和同行交流。Microsoft Ignite TheTour在大中华区有五站,分别在北京、深圳、台北、上海和香港举行。

我是微软MVP受到微软邀请,进行了3个主题的分享,分别是选择下一版 Microsoft Edge 的十大理由,面向数据的思维模式和R语言编程,基于牛顿冷却定律的热度排名算法。

目录

  1. 会议背景
  2. 我分享的主题1:选择下一版 Microsoft Edge 的十大理由
  3. 我分享的主题2:面向数据的思维模式和R语言编程
  4. 我分享的主题3:基于牛顿冷却定律的热度排名算法
  5. 会议体验和照片分享

1. 会议背景

2020年1月13-14日,Microsoft Ignite The Tour ShenZhen(MITT)深圳站吸引了近3,000名观众,Microsoft专家和技术社区为大家带来了干货满满的技术培训和讲座。本次活动的官方参会指南

深圳站13名社区讲师将带来23场演讲,微软讲师做另外97场演讲。以下是社区讲师演讲主题的分类。

本次大会我一共有3个主题:

我被照的最帅的一次。

微软公司的组合者和MVP的大合照(聚餐照)

2. 我分享的主题1:选择下一版 Microsoft Edge 的十大理由

这个主题对我来说是全新的,属于M365的产品介绍。我虽然一直在使用浏览器,但大部分时间都是被chrome占据着,并没有对Edge有太多的关注。在大会前收到了微软MVP组招募新的产品讲师,本着兴趣和好奇,接下了这个主题,PPT下载

Microsoft Edge是面向企业的一个浏览器,他的定位确实与Chrome面向个人的浏览器的定位有很大的差别。本次分享获得了一些微软官方的材料的支持,我也仔细的学习了一下,然后转换成了中文的材料。分享主要为分10个部分,也就是10大理由。

其中,第4点可信赖的服务,第5点安全性,第9点无无缝访问您所需的信息,我认为是企业浏览器独特的,区别于人个浏览器的重要的部分。Microsoft Edge的详细介绍,请查看文章选择下一版 Microsoft Edge 的十大理由。在分享后,我也已经开始尝试下一代的Edge为办公的便利性体验了。

3. 我分享的主题2:面向数据的思维模式和R语言编程

很多公司已经完成了数据的原始积累,如何让沉睡的数据发挥价值,是急需要功课的难关!

数据项目和软件项目、互联网项目都有非常大的不同,不确定性、跨学科知识点、工程落地,都是影响数据项目成功与失败的重要因素。掌握数据思维,科学的方法论,专业的团队,便利的工具,才能让数据项目走向成功。

我主要为分4个部分进行介绍:

  1. 面向数据的思维模式
  2. 如何开展一个数据项目
  3. R语言进行数据处理
  4. R语言项目案例

真实的一个数据项目,处处是坑,必经之路至少包括7个步骤:需求讨论、数据提取、数据整合、数据清洗、特征工程、模型搭建和模型评估。

这7个步骤,目前没有标准的定义,同时也会随着目标的不同,都会有一些变化,我后面会单独详细写一篇文章进行介绍,把我所经历的成功经验和失败经验都做一下总结。

3. 我分享的主题2:基于牛顿冷却定律的热度排名算法(R语言实现)

这个主题是一个比较有意思的互联网模型应用的一个主题,介绍了牛顿冷却定律在排行榜上的一种实现。我将介绍一个数据项目从研发到实施的全流程,把理论模型结合实际的场景进行应用,从理论,到公式,再到验证,然后结合实际场景,数据模拟,形成数据产品。听众可以了解到,数据项目的建立过程,思考过程,数据探索过程,帮助大家更好的使用数据。

基于牛顿冷却定律的热度排名算法(R语言实现)
我们生活中随处能看到排行榜,有图书排名,电影排名,文章排名,音乐排名,商品排名,商铺排名等等,排在TOP10名单商品,会对我们生活中的决策产生重大的影响,所以如何设计一个客观的排行榜,是一件很重要的事情。这个场景就非常时候用到算法模型,进行客观评价和排名。

牛顿冷却定律是一种热力学模型,通过温度与时间之间的函数关系,构建出了一个指数衰减的过程。物理学中热度的衰减,同样适用在排行榜的商品衰减。对基础科学的公理和定理的研究,可以极大地开阔研究数据的思路,形成跨学科的解决方案。用自然科学的基础理论,来解决跨学科的问题是一种绝佳的方法。

我主要为分4个部分进行介绍:

  1. 排名算法背景介绍
  2. 牛顿冷却定律原理
  3. 算法模拟和R语言实现
  4. 算法应用落地

最后总结,本文从场景开始,介绍了一般思路,行业思路,理论模型,数学公式,推到过程,程序实现,最后回到场景应用,形成数据研发闭环。

对基础科学的公理和定理的研究,可以极大地开阔研究数据的思路,形成跨学科的解决方案。用自然科学的基础理论,来解决跨学科的问题是一种绝佳的方法。

希望通过本次分享,给大家带入一个新的领域。

4. 会议体验和照片分享

本次大会我体会到的一些关键字:场面大(走路走断腿),微软大战略。

4.1 会议体验证和总结

深圳会展中心太大了,注册在1楼,开放式的剧院区在1楼,主题分享教室都在5楼和6楼,真是走路走断腿啊。我为找到分享的3个地方,将近花了30分钟的时间。

另外一点,一天讲三场确实是累啊!为啥把我的分享都安排到了1天了,确实是累啊!

深圳会展中心中心的全楼层,只少能容纳2万人,1层的剧院区部分!

Ai in Clond的4位MVP。

小间认真听课的观众。

MVP的小合照,可惜没有参加上。

最后,特别感谢一直在后面帮助处理各种事情的MVP助理康爽(这次没有来),下次我会记得要补上的。

MVP!Yeah!下一站就没有下一站了,这次分享后就赶上了“新型冠状病毒”,好好在家修养生息,要等疫情过去再出门了!

转载请注明出处:
http://blog.fens.me/meeting-mvp-20200114

打赏作者

R语言中文分词包jiebaR

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-word-jiebar/

jiebaR

前言

本文挖掘是数据挖掘中一个非常重要的部分,有非常广阔的使用场景,比如我们可以对新闻事件进行分析,了解国家大事;也可以对微博信息进行分析,通过社交舆情看看大家的关注点。通过文本挖掘找到文章中的隐藏信息,对文章的结构进行分析,判断是不是同一个作者写文章;同时可以对邮件分析,结合bayes算法判断哪些是垃圾邮件,哪些是有用的邮件。

本文挖掘的第一步,就是要进行分词,分词将直接影响文本挖掘的效果。R语言在分词方面有很好的支持,接下来就给大家介绍一个不错的R语言中文分词包“结巴分词”(jiebaR)。

目录

  1. jiebaR包介绍
  2. 5分钟上手
  3. 分词引擎
  4. 配置词典
  5. 停止词过滤
  6. 关键词提取

1. jiebaR包介绍

结巴分词(jiebaR),是一款高效的R语言中文分词包,底层使用的是C++,通过Rcpp进行调用很高效。结巴分词基于MIT协议,就是免费和开源的,感谢国人作者的给力支持,让R的可以方便的处理中文文本。

官方Github的地址:https://github.com/qinwf/jiebaR

本文所使用的系统环境

  • Win10 64bit
  • R: 3.2.3 x86_64-w64-mingw32/x64 b4bit

jiebaR包是在CRAN发布的标准库,安装起来非常简单,2条命令就可以了。


~ R
> install.packages("jiebaR")
> library("jiebaR")

如果想要安装开发版本,可以使用devtools来进行安装,devtools的介绍请参考文章:在巨人的肩膀前行 催化R包开发


> library(devtools)
> install_github("qinwf/jiebaRD")
> install_github("qinwf/jiebaR")
> library("jiebaR")

开发版本安装,官方建议使用Linux系统 gcc >= 4.6 编译,Windows需要安装 Rtools。

2. 5分钟上手

5分钟上手,直接看第一个例子吧,对一段文字进行分词。


> wk = worker()

> wk["我是《R的极客理想》图书作者"]
[1] "我是" "R"    "的"   "极客" "理想" "图书" "作者"

> wk["我是R语言的深度用户"]
[1] "我"   "是"   "R"    "语言" "的"   "深度" "用户"

很简单地,2行代码,就完成了中文分词。

jiebaR提供了3种分词语句的写法,例子上面的用[]符号的语法,还可以使用<=符合语法,或者使用segment()函数。虽然形式不同,但是分词效果是一样的。使用<=符号的语法,如下


> wk<='另一种符合的语法'
[1] "另"   "一种" "符合" "的"   "语法"

使用segment()函数的语法,如下


> segment( "segment()函数语句的写法" , wk )
[1] "segment" "函数"    "语句"    "的"      "写法" 

如果你觉得很神奇,想了解如何自定义操作符的,可以检查项目的源代码quick.R文件


# <= 符号定义
`<=.qseg`<-function(qseg, code){
  if(!exists("quick_worker",envir = .GlobalEnv ,inherits = F) || 
       .GlobalEnv$quick_worker$PrivateVarible$timestamp != TIMESTAMP){
    
    if(exists("qseg",envir = .GlobalEnv,inherits = FALSE ) ) 
      rm("qseg",envir = .GlobalEnv)
    
    modelpath  = file.path(find.package("jiebaR"),"model","model.rda")
    quickparam = readRDS(modelpath)
    
    if(quickparam$dict == "AUTO") quickparam$dict = DICTPATH
    if(quickparam$hmm == "AUTO") quickparam$hmm = HMMPATH
    if(quickparam$user == "AUTO") quickparam$user = USERPATH
    if(quickparam$stop_word == "AUTO") quickparam$stop_word = STOPPATH
    if(quickparam$idf == "AUTO") quickparam$idf = IDFPATH
    
    createquickworker(quickparam)
    setactive()
  } 

  //..代码省略
}

# [ 符号定义
`[.qseg`<- `<=.qseg`

我们也可以直接对文本文件进行分词,在当前目录新建一个文本文件idea.txt。


~ notepad idea.txt

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

当然,我们运行分词程序,会在当前目录生成一个新的分词结果的文件。


> wk['./idea.txt']
[1] "./idea.segment.2016-07-20_23_25_34.txt"

打开文件idea.segment.2016-07-20_23_25_34.txt,整个本文以空格进行分词。


~ notepad idea.segment.2016-07-20_23_25_34.txt

R 的 极客 理想 系列 文章 涵盖 了 R 的 思想 使用 工具 创新 等 的 一系列 要点 以 我 个人 的 学习 和 体验 去 诠释 R 的 强大 R 语言 作为 统计学 一门 语言 一直 在 小众 领域 闪耀着 光芒 直到 大 数据 的 爆发 R 语言 变成 了 一门 炙手可热 的 数据分析 的 利器 随着 越来越 多 的 工程 背景 的 人 的 加入 R 语言 的 社区 在 迅速 扩大 成长 现在 已 不仅仅 是 统计 领域 教育 银行 电商 互联网 都 在 使用 R 语言

是不是很简单,5分钟实践就能完成分词的任务。

3. 分词引擎

在调用worker()函数时,我们实际是在加载jiebaR库的分词引擎。jiebaR库提供了7种分词引擎。

  • 混合模型(MixSegment):是四个分词引擎里面分词效果较好的类,结它合使用最大概率法和隐式马尔科夫模型。
  • 最大概率法(MPSegment) :负责根据Trie树构建有向无环图和进行动态规划算法,是分词算法的核心。
  • 隐式马尔科夫模型(HMMSegment):是根据基于人民日报等语料库构建的HMM模型来进行分词,主要算法思路是根据(B,E,M,S)四个状态来代表每个字的隐藏状态。 HMM模型由dict/hmm_model.utf8提供。分词算法即viterbi算法。
  • 索引模型(QuerySegment):先使用混合模型进行切词,再对于切出来的较长的词,枚举句子中所有可能成词的情况,找出词库里存在。
  • 标记模型(tag)
  • Simhash模型(simhash)
  • 关键词模型(keywods)

如果你不太关心引擎的事,那么直接用官方推荐的混合模型(默认选择)就行了。查看worker()函数的定义。


worker(type = "mix", dict = DICTPATH, hmm = HMMPATH, user = USERPATH,
  idf = IDFPATH, stop_word = STOPPATH, write = T, qmax = 20, topn = 5,
  encoding = "UTF-8", detect = T, symbol = F, lines = 1e+05,
  output = NULL, bylines = F, user_weight = "max")

参数列表:

  • type, 引擎类型
  • dict, 系统词典
  • hmm, HMM模型路径
  • user, 用户词典
  • idf, IDF词典
  • stop_word, 关键词用停止词库
  • write, 是否将文件分词结果写入文件,默认FALSE
  • qmax, 最大成词的字符数,默认20个字符
  • topn, 关键词数,默认5个
  • encoding, 输入文件的编码,默认UTF-8
  • detect, 是否编码检查,默认TRUE
  • symbol, 是否保留符号,默认FALSE
  • lines, 每次读取文件的最大行数,用于控制读取文件的长度。大文件则会分次读取。
  • output, 输出路径
  • bylines, 按行输出
  • user_weight, 用户权重

我们在调用worker()时,就加载了分词引擎,可以打印出来,查看分词的引擎的配置。


> wk = worker()
> wk
Worker Type:  Jieba Segment

Default Method  :  mix     # 混合模型
Detect Encoding :  TRUE    # 检查编码
Default Encoding:  UTF-8   # UTF-8
Keep Symbols    :  FALSE   # 不保留符号
Output Path     :          # 输出文件目录
Write File      :  TRUE    # 写文件
By Lines        :  FALSE   # 不行输出
Max Word Length :  20      # 最大单单词长度
Max Read Lines  :  1e+05   # 最大读入文件行数

Fixed Model Components:  

$dict                      # 系统词典
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict/jieba.dict.utf8"

$user                      # 用户词典
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict/user.dict.utf8"

$hmm                       # 隐式马尔科夫模型模型
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict/hmm_model.utf8"

$stop_word                 # 停止词,无
NULL

$user_weight               # 用户词典权重
[1] "max"

$timestamp                 # 时间戳
[1] 1469027302

$default $detect $encoding $symbol $output $write $lines $bylines can be reset.

如果我们想改变分词引擎的配置项,可以在调用worker()创建分词引擎时,也可以通过wk$XX来进行设置。如果想了解wk是什么类型的对象,我们通过pryr包的otype的函数来检查wk对象的类型。关于pryr包的详细使用,请参考文章撬动R内核的高级工具包pryr


# 加载 pryr包
> library(pryr)
> otype(wk)  # 面向对象的类型检查
[1] "S3"

> class(wk)  # 查看class是属性
[1] "jiebar"  "segment" "jieba" 

4. 配置词典

对于分词的结果好坏的关键因素是词典,jiebaR默认有配置标准的词典。对于我们的使用来说,不同行业或不同的文字类型,最好用专门的分词词典。在jiebaR中通过show_dictpath()函数可以查看默认的标准词典,可以通过上一小节介绍的配置项,来指定我们自己的词典。日常对话的常用词典,比如搜狗输入法的词库。


# 查看默认的词库位置
> show_dictpath()
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict"

# 查看目录
> dir(show_dictpath())
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict"
 [1] "backup.rda"      "hmm_model.utf8"  "hmm_model.zip"  
 [4] "idf.utf8"        "idf.zip"         "jieba.dict.utf8"
 [7] "jieba.dict.zip"  "model.rda"       "README.md"      
[10] "stop_words.utf8" "user.dict.utf8" 

看到词典目录中,包括了多个文件。

  • jieba.dict.utf8, 系统词典文件,最大概率法,utf8编码的
  • hmm_model.utf8, 系统词典文件,隐式马尔科夫模型,utf8编码的
  • user.dict.utf8, 用户词典文件,utf8编码的
  • stop_words.utf8,停止词文件,utf8编码的
  • idf.utf8,IDF语料库,utf8编码的
  • jieba.dict.zip,jieba.dict.utf8的压缩包
  • hmm_model.zip,hmm_model.utf8的压缩包
  • idf.zip,idf.utf8的压缩包
  • backup.rda,无注释
  • model.rda,无注释
  • README.md,说明文件

打开系统词典文件jieba.dict.utf8,并打印前50行。


> scan(file="D:/tool/R-3.2.3/library/jiebaRD/dict/jieba.dict.utf8",
+           what=character(),nlines=50,sep='\n',
+           encoding='utf-8',fileEncoding='utf-8')
Read 50 items
 [1] "1号店 3 n"  "1號店 3 n"  "4S店 3 n"   "4s店 3 n"  
 [5] "AA制 3 n"   "AB型 3 n"   "AT&T 3 nz"  "A型 3 n"   
 [9] "A座 3 n"    "A股 3 n"    "A輪 3 n"    "A轮 3 n"   
[13] "BB机 3 n"   "BB機 3 n"   "BP机 3 n"   "BP機 3 n"  
[17] "B型 3 n"    "B座 3 n"    "B股 3 n"    "B超 3 n"   
[21] "B輪 3 n"    "B轮 3 n"    "C# 3 nz"    "C++ 3 nz"  
[25] "CALL机 3 n" "CALL機 3 n" "CD机 3 n"   "CD機 3 n"  
[29] "CD盒 3 n"   "C座 3 n"    "C盘 3 n"    "C盤 3 n"   
[33] "C語言 3 n"  "C语言 3 n"  "D座 3 n"    "D版 3 n"   
[37] "D盘 3 n"    "D盤 3 n"    "E化 3 n"    "E座 3 n"   
[41] "E盘 3 n"    "E盤 3 n"    "E通 3 n"    "F座 3 n"   
[45] "F盘 3 n"    "F盤 3 n"    "G盘 3 n"    "G盤 3 n"   
[49] "H盘 3 n"    "H盤 3 n"

我们发现系统词典每一行都有三列,并以空格分割,第一列为词项,第二列为词频,第三列为词性标记。

打开用户词典文件user.dict.utf8,并打印前50行。


> scan(file="D:/tool/R-3.2.3/library/jiebaRD/dict/user.dict.utf8",
+      what=character(),nlines=50,sep='\n',
+      encoding='utf-8',fileEncoding='utf-8')
Read 5 items
[1] "云计算"   "韩玉鉴赏" "蓝翔 nz"  "CEO"      "江大桥"  

用户词典第一行有二列,,第一列为词项,第二列为词性标记,没有词频的列。用户词典默认词频为系统词库中的最大词频。

jiebaR包关于词典词性标记,采用ictclas的标记方法。ICTCLAS 汉语词性标注集。

代码名称帮助记忆的诠释
Ag形语素形容词性语素。形容词代码为a,语素代码g前面置以A。
a形容词取英语形容词adjective的第1个字母。
ad副形词直接作状语的形容词。形容词代码a和副词代码d并在一起。
an名形词具有名词功能的形容词。形容词代码a和名词代码n并在一起。
b区别词取汉字"别"的声母。
c连词取英语连词conjunction的第1个字母。
Dg副语素副词性语素。副词代码为d,语素代码g前面置以D。
d副词取adverb的第2个字母,因其第1个字母已用于形容词。
e叹词取英语叹词exclamation的第1个字母。
f方位词取汉字"方"的声母。
g语素绝大多数语素都能作为合成词的"词根",取汉字"根"的声母。
h前接成分取英语head的第1个字母。
i成语取英语成语idiom的第1个字母。
j简称略语取汉字"简"的声母。
k后接成分
l习用语习用语尚未成为成语,有点"临时性",取"临"的声母。
m数词取英语numeral的第3个字母,n,u已有他用。
Ng名语素名词性语素。名词代码为n,语素代码g前面置以N。
n名词取英语名词noun的第1个字母。
nr人名名词代码n和"人(ren)"的声母并在一起。
ns地名名词代码n和处所词代码s并在一起。
nt机构团体"团"的声母为t,名词代码n和t并在一起。
nz其他专名"专"的声母的第1个字母为z,名词代码n和z并在一起。
o拟声词取英语拟声词onomatopoeia的第1个字母。
p介词取英语介词prepositional的第1个字母。
q量词取英语quantity的第1个字母。
r代词取英语代词pronoun的第2个字母,因p已用于介词。
s处所词取英语space的第1个字母。
Tg时语素时间词性语素。时间词代码为t,在语素的代码g前面置以T。
t时间词取英语time的第1个字母。
u助词取英语助词auxiliary 的第2个字母,因a已用于形容词。
Vg动语素动词性语素。动词代码为v。在语素的代码g前面置以V。
v动词取英语动词verb的第一个字母。
vd副动词直接作状语的动词。动词和副词的代码并在一起。
vn名动词指具有名词功能的动词。动词和名词的代码并在一起。
w标点符号
x非语素字非语素字只是一个符号,字母x通常用于代表未知数、符号。
y语气词取汉字"语"的声母。
z状态词取汉字"状"的声母的前一个字母。

下面我们自定义一个用户词典,来试试效果。编写词典文件,user.utf8。


~ notepad user.utf8

R语言
R的极客理想
大数据
数据

使用我们的自定义的用户词典,对刚才的文本再进行分词。


> wk = worker(user='user.utf8')
> wk['./idea.txt']
[1] "./idea.segment.2016-07-21_11_14_24.txt"

对比2次产生的分词结果,idea.segment.2016-07-20_23_25_34.txt 和 idea.segment.2016-07-21_11_14_24.txt。

jiebaR-cut

在实际使用中,jiebaR默认提供的用户词典只有5个单词,太简单了,肯定是不够用的。我们可以用搜狗词典,来丰富用户自己的词库。接下来,让我们配置搜狗词典。你需要安装一个搜狗输入法,具体的安装过程不再解释。

我安装的是搜狗五笔输入法,找到搜狗的安装目录,并找到词典文件。我的搜狗词典,在下面的安装位置。


C:\Program Files (x86)\SogouWBInput\2.1.0.1288\scd\17960.scel

把17960.scel文件复制到自己的项目目录里,用文本编辑器打开文件,发现是二进制的。那么我需要用工具进行转换,把二进制的词典转成我们可以使用的文本文件。jiebaR包的作者,同时开发了一个cidian项目,可以转换搜狗的词典,那么我们只需要安装cidian包即可。

安装cidian项目


> install.packages("devtools")
> install.packages("stringi")
> install.packages("pbapply")
> install.packages("Rcpp")
> install.packages("RcppProgress")
> library(devtools)
> install_github("qinwf/cidian")
> library(cidian)

转换二进制词典到文本文件。


# 转换
> decode_scel(scel = "./17960.scel",cpp = TRUE)
output file: ./17960.scel_2016-07-21_00_22_11.dict

# 查看生成的词典文件
> scan(file="./17960.scel_2016-07-21_00_22_11.dict",
+      what=character(),nlines=50,sep='\n',
+      encoding='utf-8',fileEncoding='utf-8')
Read 50 items
 [1] "阿坝州 n"         "阿百川 n"         "阿班 n"          
 [4] "阿宾 n"           "阿波菲斯 n"       "阿不都热希提 n"  
 [7] "阿不都西库尔 n"   "阿不力克木 n"     "阿尔姆格伦 n"    
[10] "阿尔沙文 n"       "阿肥星 n"         "阿菲正传 n"      
[13] "阿密特 n"         "阿穆 n"           "阿穆隆 n"        
[16] "阿帕鲁萨镇 n"     "阿披实 n"         "阿衰 n"          
[19] "阿霞 n"           "艾奥瓦 n"         "爱不疚 n"        
[22] "爱的错位 n"       "爱得得体 n"       "爱的火焰 n"      
[25] "爱的流刑地 n"     "爱得起 n"         "埃夫隆 n"        
[28] "爱搞网 n"         "爱国红心 n"       "爱呼 n"          
[31] "爱就宅一起 n"     "埃克希儿 n"       "爱没有错 n"      
[34] "埃蒙斯 n"         "爱奴新传 n"       "爱起点 n"        
[37] "爱情的牙齿 n"     "爱情海滨 n"       "爱情节 n"        
[40] "爱情美的样子 n"   "爱情无限谱 n"     "爱情占线 n"      
[43] "爱情转移 n"       "爱情左灯右行 n"   "爱上你是一个错 n"
[46] "矮哨兵 n"         "爱是妥协 n"       "爱似水仙 n"      
[49] "爱太痛 n"         "爱无界 n"    

接下来,直接把搜狗词典配置到我们的分词库中,就可以直接使用了。把搜狗词典文件改名,从17960.scel_2016-07-21_00_22_11.dict到user.dict.utf8,然后替换D:\tool\R-3.2.3\library\jiebaRD\dict目录下面的user.dict.utf8。这样默认的用户词典,就是搜狗词典了。很酷吧!

5. 停止词过滤

停止词就是分词过程中,我们不需要作为结果的词,像英文的语句中有很多的a,the,or,and等,中文语言中也有很多,比如 的,地,得,我,你,他。这些词因为使用频率过高,会大量出现在一段文本中,对于分词后的结果,在统计词频的时候会增加很多的噪音,所以我们通常都会将这些词进行过滤。

在jiebaR中,过滤停止词有2种方法,一种是通过配置stop_word文件,另一种是使用filter_segment()函数。

首先我们先来看,通过配置stop_word文件的方法。新建一个stop_word.txt文件。


~ notepad stop_word.txt

我
我是

加载分词引擎,并配置停止词过滤。


> wk = worker(stop_word='stop_word.txt')
> segment<-wk["我是《R的极客理想》图书作者"]
> segment
[1] "R"    "的"   "极客" "理想" "图书" "作者"

上面的文本,我们把"我是"通过停止词进行了过滤。如果还想过滤“作者”一词,可以动态的调用filter_segment()函数。


> filter<-c("作者")
> filter_segment(segment,filter)
[1] "R"    "的"   "极客" "理想" "图书"

6. 关键词提取

关键词提取是文本处理非常重要的一个环节,一个经典算法是TF-IDF算法。其中,TF(Term Frequency)代表词频,IDF(Inverse Document Frequency)表示逆文档频率。如果某个词在文章中多次出现,而且不是停止词,那么它很可能就反应了这段文章的特性,这就是我们要找的关键词。再通过IDF来算出每个词的权重,不常见的词出现的频率越高,则权重越大。计算TF-IDF的公式为:

TF-IDF = TF(词频) * 逆文档频率(IDF)

对文档中每个词计算TF-IDF的值,把结果从大到小排序,就得到了这篇文档的关键性排序列表。关于IF-IDF的解释,参考了文章TF-IDF与余弦相似性的应用(一):自动提取关键词

jiebaR包的关键词提取提取的实现,也是使用了TF-IDF的算法。在安装目录中的idf.utf8文件,为IDF的语料库。查看idf.utf8内容。


> scan(file="D:/tool/R-3.2.3/library/jiebaRD/dict/idf.utf8",
+      what=character(),nlines=50,sep='\n',
+      encoding='utf-8',fileEncoding='utf-8')
Read 50 items
 [1] "劳动防护 13.900677652"      "生化学 13.900677652"       
 [3] "奥萨贝尔 13.900677652"      "考察队员 13.900677652"     
 [5] "岗上 11.5027823792"         "倒车档 12.2912397395"      
 [7] "编译 9.21854642485"         "蝶泳 11.1926274509"        
 [9] "外委 11.8212361103"         "故作高深 11.9547675029"    
[11] "尉遂成 13.2075304714"       "心源性 11.1926274509"      
[13] "现役军人 10.642581114"      "杜勃留 13.2075304714"      
[15] "包天笑 13.900677652"        "贾政陪 13.2075304714"      
[17] "托尔湾 13.900677652"        "多瓦 12.5143832909"        
[19] "多瓣 13.900677652"          "巴斯特尔 11.598092559"     
[21] "刘皇帝 12.8020653633"       "亚历山德罗夫 13.2075304714"
[23] "社会公众 8.90346537821"     "五百份 12.8020653633"      
[25] "两点阈 12.5143832909"       "多瓶 13.900677652"         
[27] "冰天 12.2912397395"         "库布齐 11.598092559"       
[29] "龙川县 12.8020653633"       "银燕 11.9547675029"        
[31] "历史风貌 11.8212361103"     "信仰主义 13.2075304714"    
[33] "好色 10.0088573539"         "款款而行 12.5143832909"    
[35] "凳子 8.36728816325"         "二部 9.93038573842"        
[37] "卢巴 12.1089181827"         "五百五 13.2075304714"      
[39] "畅叙 11.598092559"          "吴栅子 13.2075304714"      
[41] "智力竞赛 13.900677652"      "库邦 13.2075304714"        
[43] "非正义 11.3357282945"       "编订 10.2897597393"        
[45] "悲号 12.8020653633"         "陈庄搭 13.2075304714"      
[47] "二郎 9.62401153296"         "电光石火 11.8212361103"    
[49] "抢球 11.9547675029"         "南澳大利亚 10.9562386728"  

idf.utf8文件每一行有2列,第一列是词项,第二列为权重。然后,我通过计算文档的词频(TF),与语料库的IDF值相乘,就可以得到TF-IDF值,从而提取文档的关键词。

比如,我们对下面的文本内容进行关键词的提取。


> wk = worker()
> segment<-wk["R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。"]

# 计算词频
> freq(segment)
     char freq
1    创新    1
2      了    1
3    文章    1
4    强大    1
5       R    3
6    个人    1
7      的    5
8    诠释    1
9      和    1
10 一系列    1
11   使用    1
12     以    1
13     等    1
14   极客    1
15   理想    1
16   思想    1
17   涵盖    1
18   系列    1
19     去    1
20     我    1
21   工具    1
22   学习    1
23   体验    1
24   要点    1

# 取TF-IDF的前5的关键词
> keys = worker("keywords",topn=5)

# 计算关键词
> vector_keywords(segment,keys)
11.7392 8.97342 8.23425  8.2137 7.43298 
 "极客"  "诠释"  "要点"  "涵盖"  "体验" 

使用jiebaR包处理分词确实简单,几行的代码就能实现分词的各种算法操作。有了这个工具,我们就可以文档中,发现各种语言规则进行文本挖掘了。下篇文章让我们挖掘一下上市公司的公告吧,说不定能发现什么市场规则。

本文只是抛砖引玉地介绍了jiebaR包的使用方法,详细使用操作,请参考包作者的官方介绍。再次感谢jiebaR作者@qinwenfeng,为R语言在中文分词中提供了一套非常不错的工具包!

转载请注明出处:
http://blog.fens.me/r-word-jiebar/

打赏作者

2016天善智能交流会第22场: R语言为量化而生

跨界知识聚会系列文章,“知识是用来分享和传承的”,各种会议、论坛、沙龙都是分享知识的绝佳场所。我也有幸作为演讲嘉宾参加了一些国内的大型会议,向大家展示我所做的一些成果。从听众到演讲感觉是不一样的,把知识分享出来,你才能收获更多。

关于作者

  • 张丹(Conan), 程序员Java,R,Nodejs
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/meeting-hellobi-20160701/

meeting-hellobi

前言

感谢天善智能社区的邀请,有幸参加每周一期的跟数据有关的行业、工具、技术的交流盛宴,活动的口号是“Friday BI Fly 周五BI飞起来”。

目录

  1. 我的分享主题:R语言为量化而生
  2. 会议体验
  3. 自由讨论

1. 我的分享主题:R语言为量化而生

本次分享的主题 R语言为量化而生,主要内容来自我的一篇博客文章:R语言为量化而生。希望能够解释清楚,在量化投资中为什么要用R语言。从程序员的角度看,C++,Java,Python, C#都是可行方案;从数据人员的角度看,Excel, SAS, Matlab更是不错的。那么为什么是R语言呢,R语言的优势在哪里体现?

这类的问题,总是会被问到。那么答案,就在于你对量化这件事情的了解,和对各种编程语言的理解。最近3年,互联网在量化领域的大发展,以Quantopian为代表的在线策略研发平台,用Python做为核心语言,国内同样支持Python的平台也有 优矿聚宽米筐。这些平台主是面向程序员群体的平台,希望通过挖掘草根明星,来推动量化的发展。传统的量化交易软件,像文华MC, TB, TS 都有自己一套的脚本化的编程语言。有实力的专业团队,通常会自成体系的独立开发一套自己的系统。如果面向更广泛的人群,最常用的方法就是Wind导数据,Excel中拉个表出来。

所以,其实用什么语言不重要,关键是怎么理解做量化这件事情。那么R语言的天生优势就是数学计算,数据处理,免费开源,大量支持库。试试吧,你一定会喜欢的。

2. 会议体验

本次分享受天善智能社区的邀请,我真的非常高兴。天善智能是新一代的商业智能和大数据的垂直社区,聚集了大量的数据分析从业人员。活动介绍,https://ask.hellobi.com/blog/tianshansoft/4229。 本次活动同时有30个微信群进行直播,参加的人员,至少有2000人以上。可以天善智能社区,在行业的影响力是非常大的。

发个截图,体会一下微信同步直播的震撼吧!

wx

本此的分享基于微信的直播,我也第一次体验,要用纯文字的方式来进行介绍。想把一个事情说清楚,又增加了不少的难度。由于不能分享屏幕,代码部分会通过图片截屏。

本次活动的总结,https://ask.hellobi.com/blog/tianshansoft/4271,感谢天善社区的工作人员进行整理。

远程分享,就是没能与大家合照,有点遗憾!!贴张自己的照片吧。

01

3. 自由讨论

分享后,很多朋友都对于R语言都是非常的好奇,提了很多的问题,用户的参与性非常强。下列直接贴出用户的问题和我的回复。

1、替新手问一个,请教一下,R语言的数据分析应该从哪方面入手练习啊?因为目前工作上不是用R的,看完书之后想具体去试一下。

张丹: R其实上手很快,找一本书,认真操作练习一遍就上手了。

2、玉琴:不建议用for loop的原因是考虑到性能问题吗

张丹:for loop是调用的R的循环库,apply是调用C的循环库,性能差距还是很大的

3、来自20群的提问:提个问题,微软对R的收购会对R语言的发展产生什么影响?

张丹:我觉得这是正向发展的,是好事情。大公司看到了R的潜力!

4、尚林栋:R语言金融建模的具体步骤能说一下吗

金融建模的具体步骤,你可以参考这篇文章,http://blog.fens.me/finance-stock-ma/

5、刘嘉丰Alan:丹哥,现在有很多量化平台,提供打包好的函数,在线回测,和自己造轮子拿R语言相比,您觉得各有什么优势呢?

张丹:R的优势就是在数学计算,数据处理上。行业标准还没有统一,所以不一定在线平台的轮子就一定好用。但另外,我们从开发或使用的角度,更多的用到的R包,都是RStudio公司的产品,我觉得是RStudio在推动R的整个的进化过程。

6、我也觉得r语言不错,但经常想不到商业场景,到现在,我只是用它统计考勤,各种绩效kpi,每月算一次奖金,已经这样过去2年了,r语言路在何方哪?

张丹: 你所说的统计,只能说简单计数。比如,你要预测下个月的考勤情况,从而设计预算方案。你可能就需要做个回归分析,这时R就能给你很大的帮助了。生活和工作中,随处都是数据分析的场景。

7、Allen:r在拟合上感觉比python用起来更爽一些,其返回的结果较多

张丹:那么R和python比,R更面向数据,特别是对于没有编程基础的人。PYTHON,还是程序语言,还要了解程序结构,程序架构,代码量不会少。

有IT背景程序员,可能更倾向于PYTHON;如果没有IT背景,R更容易上手。

8、越中女儿:请教一个问题:quantmod对美股的实时接口很好用,对A股不支持,且A股基本面数据才更新到2013.09,请问有好用的ETL包么,类似于python的tushare那样对A股友好的,各种etl啊清洗的脏活累活感觉python更好啊,R就是安安静静做做统计,玩玩图形。

张丹: quantmod使用的是yahoo等国外的数据源,这些数据源本身没有A股数据,如果需要A股数据,用tushare还是不错的。 R特有的data.frame,matrix 等类型和操作方法,在python也需要单独去实现。

9、柠檬味的香草:最近想研究一些互联网文本数据与指数或各股走势的关系,但是在使用R语言处理文本数据不是很方便,丹哥可有一些强大的library推荐,对于非结构,文本数据的处理。

张丹:“尽量使用向量计算或矩阵计算的计算方法”,可以这样理解,对于一个二维结构,for需要2次,0(N^2)的时间复杂度。如果我们把数据,直接就按矩阵存储, 你让矩阵里的每个点都加1, 只需要算一次。Hadley提供的包,源代码我都看过,写很棒,也很实用。

r在拟合上感觉比python用起来更爽一些,其返回的结果较多

其实R有很多的第三方的包,已经有了大量的算法包,而其他语言相对较少。只是我们平时接触的不多,所以觉得用不到。R有大量的统计包,你可以从官方网站找到,输出的结果,大部分也都是统计的结果。

R所支持的行业领域,非常广泛。而工程的语言,不会做细粒度的区分,只是通用的解决方法。

10、郑州—金融数据:python有pandas.DataFrame,pandas应该是第三方的数据库结构吧?R的data.frame是内置的。

张丹:pandas.DataFrame,在底层处理,还需要对原PYTHON的数据结构做映射。当然他可以解决的很好,但你看到的内存结构,可能并不是真正的内存结构。

R内置数据类型,就可以理解是内存结构。不需要再考虑转换了。找一个自己熟悉的语言,大多数的功能,每种语言都是能实现。只有很细的领域,才会进一步区分。

11、RHaoop采用分布式并行计算,那请问如何解决需要嵌套循环的算法。

张丹:对于基于hadoop大数据的MR计算,建议做数学变成,通过数学的角度处理。我写过2个例子,一个是pagerank, 一个是itemcf。

12、@柠檬味的香草:想听听丹哥对传统数据挖掘转量化投资的建议。比如前景?竞争力?

张丹:量化投资,其实是IT人都想转的行业。你写的代码,不是通过工资来赚钱,而直接通过交易赚钱,代码的效用是最大化的。但这个行业竞争很大,聪明人都在这里,要么你的技术牛,要么你了解市场,要么的算法是独特的,不然也很难。

JhT: 做量化交易和策略的都是高智商的

越中女儿:我觉得量化对金融市场的理解比对技术本身更重要,R的需求应该会很快凸显出来。因为数据基础都有了,后面就是差会分析的人了。通常懂数据分析的程序员,比纯程序员待遇高。

13、老师,有好的spark或者hadoop入门的书吗,计算机能力弱和java不懂啊

张丹:hadoop有很多书了,我当初看的是 权威指南。spark的书不了解,我的是网上文档。

14、@Mia.W 学RHadoop需要对Hadoop或Mapreduce了解到什么程度,需要从头学hadoop或java吗

张丹:hadoop的MR的原理要了解,找到懂JAVA的同事,帮你把环境搭好。

15、@JhT 我是刚进来的,R的优势是什么?

张丹:R是免费开源的,CRAN上有8000多个包,遍布各行各业。R语言的3个特性,数学计算,数据建模,可视化。

16、@郑州—金融数据个人感觉商业上matlab比R和python支持度都要好,不管是分析,统计,挖掘还是量化方便,收费的毕竟是收费的

张丹:有商业推动,当然要比免费的好了。不过,像SAS和Matlab也在打通和R的接口,毕竟由全球第三方贡献包,要比一家公司提供的包要多很多的。

17、@越中女儿 有用R做过实盘风控么

张丹:有做,其实不太复杂。你把需要的实时数据,都同步存到redis中,用R在秒级调reids取数据,计算完成再写回去。

18、@Jason.k计算机8g内存,数据虽然行数不多,但是很多列,所以数据csv格式大小会高达几个G,这个规模数据量,内存应该是不够的。

张丹:R的机制,会把数据一次性加载到内存中。就算能读到内存,每次计算时,也会有中间变量,所以你的基础内存是不够的。而且对于win性能会更差。

最后,再次感谢 天善社区的小伙伴们的努力,谢谢大家!

转载请注明出处:
http://blog.fens.me/meeting-hellobi-20160701/

打赏作者

超高性能数据处理包data.table

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-data-table/

datatable-title

前言

在R语言中,我们最常用的数据类型是data.frame,绝大多数的数据处理的操作都是围绕着data.frame结构来做的。用data.frame可以很方便的进行数据存储和数据查询,配合apply族函数对数据循环计算,也可也用plyr, reshape2, melt等包对数据实现切分、分组、聚合等的操作。在数据量不太大的时候,使用起来很方便。但是,用data.frame结构处理数据时并不是很高效,特别是在稍大一点数据规模的时候,就会明显变慢。

data.table其实提供了一套和data.frame类似的功能,特别增加了索引的设置,让数据操作非常高效,可能会提升1-2数量级。本章就将data.table包的使用方法。

目录

  1. data.table包介绍
  2. data.table包的使用
  3. data.table包性能对比

1. data.table包介绍

data.table包是一个data.frame的扩展工具集,可以通过自定义keys来设置索引,实现高效的数据索引查询、快速分组、快速连接、快速赋值等数据操作。data.table主要通过二元检索法大大提高数据操作的效率,它也兼容适用于data.frame的向量检索法。同时,data.table对于大数据的快速聚合也有很好的效果,官方介绍说对于 100GB规模内存数据处理,运行效率还是很好的。那么,就让我们试验一下吧。

data.table项目地址:https://cran.r-project.org/web/packages/data.table/

本文所使用的系统环境

  • Win10 64bit
  • R: 3.2.3 x86_64-w64-mingw32/x64 b4bit

data.table包是在CRAN发布的标准库,安装起来非常简单,2条命令就可以了。


~ R
> install.packages("data.table")
> library(data.table)

2. data.table包的使用

接下来,开始用data.table包,并熟悉一下data.table包的基本操作。

2.1 用data.table创建数据集

通常情况,我们用data.frame创建一个数据集时,可以使用下面的语法。


# 创建一个data.frame数据框
> df<-data.frame(a=c('A','B','C','A','A','B'),b=rnorm(6))
> df
  a          b
1 A  1.3847248
2 B  0.6387315
3 C -1.8126626
4 A -0.0265709
5 A -0.3292935
6 B -1.0891958

对于data.table来说,创建一个数据集是和data.frame同样语法。


# 创建一个data.table对象
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a           b
1: A  0.09174236
2: B -0.84029180
3: C -0.08157873
4: A -0.39992084
5: A -1.66034154
6: B -0.33526447

检查df, dt两个对象的类型,可以看到data.table是对data.frame的扩展类型。


# data.frame类型
> class(df)
[1] "data.frame"

# data.table类型
> class(dt)
[1] "data.table" "data.frame"

如果data.table仅仅是对data.frame的做了S3的扩展类型,那么data.table是不可能做到对data.frame从效率有极大的改进的。为了验证,我们需要检查一下data.table代码的结构定义。


# 打印data.table函数定义
> data.table
function (..., keep.rownames = FALSE, check.names = FALSE, key = NULL) 
{
    x <- list(...)
    if (!.R.listCopiesNamed) 
        .Call(CcopyNamedInList, x)
    if (identical(x, list(NULL)) || identical(x, list(list())) || 
        identical(x, list(data.frame(NULL))) || identical(x, 
        list(data.table(NULL)))) 
        return(null.data.table())
    tt <- as.list(substitute(list(...)))[-1L]
    vnames = names(tt)
    if (is.null(vnames)) 
        vnames = rep.int("", length(x))
    vnames[is.na(vnames)] = ""
    novname = vnames == ""
    if (any(!novname)) {
        if (any(vnames[!novname] == ".SD")) 
            stop("A column may not be called .SD. That has special meaning.")
    }
    for (i in which(novname)) {
        if (is.null(ncol(x[[i]]))) {
            if ((tmp <- deparse(tt[[i]])[1]) == make.names(tmp)) 
                vnames[i] <- tmp
        }
    }
    tt = vnames == ""
    if (any(tt)) 
        vnames[tt] = paste("V", which(tt), sep = "")
    n <- length(x)
    if (n < 1L) 
        return(null.data.table())
    if (length(vnames) != n) 
        stop("logical error in vnames")
    vnames <- as.list.default(vnames)
    nrows = integer(n)
    numcols = integer(n)
    for (i in seq_len(n)) {
        xi = x[[i]]
        if (is.null(xi)) 
            stop("column or argument ", i, " is NULL")
        if ("POSIXlt" %chin% class(xi)) {
            warning("POSIXlt column type detected and converted to POSIXct. We do not recommend use of POSIXlt at all because it uses 40 bytes to store one date.")
            x[[i]] = as.POSIXct(xi)
        }
        else if (is.matrix(xi) || is.data.frame(xi)) {
            xi = as.data.table(xi, keep.rownames = keep.rownames)
            x[[i]] = xi
            numcols[i] = length(xi)
        }
        else if (is.table(xi)) {
            x[[i]] = xi = as.data.table.table(xi, keep.rownames = keep.rownames)
            numcols[i] = length(xi)
        }
        nrows[i] <- NROW(xi)
        if (numcols[i] > 0L) {
            namesi <- names(xi)
            if (length(namesi) == 0L) 
                namesi = rep.int("", ncol(xi))
            namesi[is.na(namesi)] = ""
            tt = namesi == ""
            if (any(tt)) 
                namesi[tt] = paste("V", which(tt), sep = "")
            if (novname[i]) 
                vnames[[i]] = namesi
            else vnames[[i]] = paste(vnames[[i]], namesi, sep = ".")
        }
    }
    nr <- max(nrows)
    ckey = NULL
    recycledkey = FALSE
    for (i in seq_len(n)) {
        xi = x[[i]]
        if (is.data.table(xi) && haskey(xi)) {
            if (nrows[i] < nr) 
                recycledkey = TRUE
            else ckey = c(ckey, key(xi))
        }
    }
    for (i in which(nrows < nr)) {
        xi <- x[[i]]
        if (identical(xi, list())) {
            x[[i]] = vector("list", nr)
            next
        }
        if (nrows[i] == 0L) 
            stop("Item ", i, " has no length. Provide at least one item (such as NA, NA_integer_ etc) to be repeated to match the ", 
                nr, " rows in the longest column. Or, all columns can be 0 length, for insert()ing rows into.")
        if (nr%%nrows[i] != 0L) 
            warning("Item ", i, " is of size ", nrows[i], " but maximum size is ", 
                nr, " (recycled leaving remainder of ", nr%%nrows[i], 
                " items)")
        if (is.data.frame(xi)) {
            ..i = rep(seq_len(nrow(xi)), length.out = nr)
            x[[i]] = xi[..i, , drop = FALSE]
            next
        }
        if (is.atomic(xi) || is.list(xi)) {
            x[[i]] = rep(xi, length.out = nr)
            next
        }
        stop("problem recycling column ", i, ", try a simpler type")
        stop("argument ", i, " (nrow ", nrows[i], ") cannot be recycled without remainder to match longest nrow (", 
            nr, ")")
    }
    if (any(numcols > 0L)) {
        value = vector("list", sum(pmax(numcols, 1L)))
        k = 1L
        for (i in seq_len(n)) {
            if (is.list(x[[i]]) && !is.ff(x[[i]])) {
                for (j in seq_len(length(x[[i]]))) {
                  value[[k]] = x[[i]][[j]]
                  k = k + 1L
                }
            }
            else {
                value[[k]] = x[[i]]
                k = k + 1L
            }
        }
    }
    else {
        value = x
    }
    vnames <- unlist(vnames)
    if (check.names) 
        vnames <- make.names(vnames, unique = TRUE)
    setattr(value, "names", vnames)
    setattr(value, "row.names", .set_row_names(nr))
    setattr(value, "class", c("data.table", "data.frame"))
    if (!is.null(key)) {
        if (!is.character(key)) 
            stop("key argument of data.table() must be character")
        if (length(key) == 1L) {
            key = strsplit(key, split = ",")[[1L]]
        }
        setkeyv(value, key)
    }
    else {
        if (length(ckey) && !recycledkey && !any(duplicated(ckey)) && 
            all(ckey %in% names(value)) && !any(duplicated(names(value)[names(value) %in% 
            ckey]))) 
            setattr(value, "sorted", ckey)
    }
    alloc.col(value)
}
<bytecode: 0x0000000017bfb990>
<environment: namespace:data.table>

从上面的整个大段代码来看,data.table的代码定义中并没有使用data.frame结构的依赖的代码,data.table都在自己函数定义中做的数据处理,所以我们可以确认data.table和data.frame的底层结果是不一样的。

那么为什么从刚刚用class函数检查data.table对象时,会看到data.table和data.frame的扩展关系呢?这里就要了解R语言中对于S3面向对象系统的结构设计了,关于S3的面向对象设计,请参考文章R语言基于S3的面向对象编程

从上面代码中,倒数第17行找到 setattr(value, "class", c("data.table", "data.frame")) 这行,发现这个扩展的定义是作者主动设计的,那么其实就可以理解为,data.table包的作者希望data.table使用起来更像data.frame,所以通过一些包装让使用者无切换成本的。

2.2 data.table和data.frame相互转换

如果想把data.frame对象和data.table对象进行转换,转换的代码是非常容易的,直接转换就可以了。

从一个data.frame对象转型到data.table对象。


# 创建一个data.frame对象
> df<-data.frame(a=c('A','B','C','A','A','B'),b=rnorm(6))

# 检查类型
> class(df)
[1] "data.frame"

# 转型为data.table对象
> df2<-data.table(df)

# 检查类型
> class(df2)
[1] "data.table" "data.frame"

从一个data.table对象转型到data.frame对象。


# 创建一个data.table对象
> dt <- data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))

# 检查类型
> class(dt)
[1] "data.table" "data.frame"

# 转型为data.frame对象
> dt2<-data.frame(dt)

# 检查类型
> class(dt2)
[1] "data.frame"

2.3 用data.table进行查询

由于data.table对用户使用上是希望和data.frame的操作尽量相似,所以适用于data.frame的查询方法基本都适用于data.table,同时data.table自己具有的一些特性,提供了自定义keys来进行高效的查询。

下面先看一下,data.table基本的数据查义方法。


# 创建一个data.table对象
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a          b
1: A  0.7792728
2: B  1.4870693
3: C  0.9890549
4: A -0.2769280
5: A -1.3009561
6: B  1.1076424

按行或按列查询


# 取第二行的数据
> dt[2,]
   a        b
1: B 1.487069

# 不加,也可以
> dt[2]
   a        b
1: B 1.487069


# 取a列的值
> dt$a
[1] "A" "B" "C" "A" "A" "B"

# 取a列中值为B的行
> dt[a=="B",]
   a        b
1: B 1.487069
2: B 1.107642

# 取a列中值为B的行的判断
> dt[,a=='B']
[1] FALSE  TRUE FALSE FALSE FALSE  TRUE

# 取a列中值为B的行的索引
> which(dt[,a=='B'])
[1] 2 6

上面的操作,不管是用索引值,== 和 $ 都是data.frame操作一样的。下面我们取data.table特殊设计的keys来查询。


# 设置a列为索引列
> setkey(dt,a)

# 打印dt对象,发现数据已经按照a列字母对应ASCII码值进行了排序。
> dt
   a          b
1: A  0.7792728
2: A -0.2769280
3: A -1.3009561
4: B  1.4870693
5: B  1.1076424
6: C  0.9890549

按照自定义的索引进行查询。


# 取a列中值为B的行
> dt["B",]
   a        b
1: B 1.487069
2: B 1.107642

# 取a列中值为B的行,并保留第一行
> dt["B",mult="first"]
   a        b
1: B 1.487069

# 取a列中值为B的行,并保留最后一行
> dt["B",mult="last"]
   a        b
1: B 1.107642

# 取a列中值为b的行,没有数据则为NA
> dt["b"]
   a  b
1: b NA

从上面的代码测试中我们可以看出,在定义了keys后,我们要查询的时候就不用再指定列了,默认会把方括号中的第一位置留给keys,作为索引匹配的查询条件。从代码的角度,又节省了一个变量定义的代码。同时,可以用mult参数,对数据集增加过滤条件,让代码本身也变得更高效。如果查询的值,不是索引列包括的值,则返回NA。

2.4 对data.table对象进行增、删、改操作

给data.table对象增加一列,可以使用这样的格式 data.table[, colname := var1]。


# 创建data.table对象
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a           b
1: A  1.51765578
2: B  0.01182553
3: C  0.71768667
4: A  0.64578235
5: A -0.04210508
6: B  0.29767383

# 增加1列,列名为c
> dt[,c:=b+2]
> dt
   a           b        c
1: A  1.51765578 3.517656
2: B  0.01182553 2.011826
3: C  0.71768667 2.717687
4: A  0.64578235 2.645782
5: A -0.04210508 1.957895
6: B  0.29767383 2.297674

# 增加2列,列名为c1,c2
> dt[,`:=`(c1 = 1:6, c2 = 2:7)]
> dt
   a          b        c c1 c2
1: A  0.7545555 2.754555  1  2
2: B  0.5556030 2.555603  2  3
3: C -0.1080962 1.891904  3  4
4: A  0.3983576 2.398358  4  5
5: A -0.9141015 1.085899  5  6
6: B -0.8577402 1.142260  6  7

# 增加2列,第2种写法
> dt[,c('d1','d2'):=list(1:6,2:7)]
> dt
   a          b        c c1 c2 d1 d2
1: A  0.7545555 2.754555  1  2  1  2
2: B  0.5556030 2.555603  2  3  2  3
3: C -0.1080962 1.891904  3  4  3  4
4: A  0.3983576 2.398358  4  5  4  5
5: A -0.9141015 1.085899  5  6  5  6
6: B -0.8577402 1.142260  6  7  6  7

给data.table对象删除一列时,就是给这列赋值为空,使用这样的格式 data.table[, colname := NULL]。我们继续使用刚才创建的dt对象。


# 删除c1列
> dt[,c1:=NULL]
> dt
   a          b        c c2 d1 d2
1: A  0.7545555 2.754555  2  1  2
2: B  0.5556030 2.555603  3  2  3
3: C -0.1080962 1.891904  4  3  4
4: A  0.3983576 2.398358  5  4  5
5: A -0.9141015 1.085899  6  5  6
6: B -0.8577402 1.142260  7  6  7

# 同时删除d1,d2列
> dt[,c('d1','d2'):=NULL]
> dt
   a          b        c c2
1: A  0.7545555 2.754555  2
2: B  0.5556030 2.555603  3
3: C -0.1080962 1.891904  4
4: A  0.3983576 2.398358  5
5: A -0.9141015 1.085899  6
6: B -0.8577402 1.142260  7

修改data.table对象的值,就是通过索引定位后进行值的替换,通过这样的格式 data.table[condition, colname := 0]。我们继续使用刚才创建的dt对象。


# 给b赋值为30
> dt[,b:=30]
> dt
   a  b        c c2
1: A 30 2.754555  2
2: B 30 2.555603  3
3: C 30 1.891904  4
4: A 30 2.398358  5
5: A 30 1.085899  6
6: B 30 1.142260  7

# 对a列值为B的行,c2列值值大于3的行,的b列赋值为100
> dt[a=='B' & c2>3, b:=100]
> dt
   a   b        c c2
1: A  30 2.754555  2
2: B  30 2.555603  3
3: C  30 1.891904  4
4: A  30 2.398358  5
5: A  30 1.085899  6
6: B 100 1.142260  7

# 还有另一种写法
> dt[,b:=ifelse(a=='B' & c2>3,50,b)]
> dt
   a  b        c c2
1: A 30 2.754555  2
2: B 30 2.555603  3
3: C 30 1.891904  4
4: A 30 2.398358  5
5: A 30 1.085899  6
6: B 50 1.142260  7

2.5 data.table的分组计算

基于data.frame对象做分组计算时,要么使用apply函数自己处理,要么用plyr包的分组计算功能。对于data.table包本身就支持了分组计算,很像SQL的group by这样的功能,这是data.table包主打的优势。

比如,按a列分组,并对b列按分组求和。


# 创建数据
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a          b
1: A  1.4781041
2: B  1.4135736
3: C -0.6593834
4: A -0.1231766
5: A -1.7351749
6: B -0.2528973

# 对整个b列数据求和
> dt[,sum(b)]
[1] 0.1210455

# 按a列分组,并对b列按分组求和
> dt[,sum(b),by=a]
   a         V1
1: A -0.3802474
2: B  1.1606763
3: C -0.6593834

2.6 多个data.table的连接操作

在操作数据的时候,经常会出现2个或多个数据集通过一个索引键进行关联,而我们的算法要把多种数据合并到一起再进行处理,那么这个时候就会用的数据的连接操作,类似关系型数据库的左连接(LEFT JOIN)。

举个例子,学生考试的场景。按照ER设计方法,我们通常会按照实体进行数据划分。这里存在2个实体,一个是学生,一个是成绩。学生实体会包括,学生姓名等的基本资料,而成绩实体会包括,考试的科目,考试的成绩。

假设有6个学生,分别参加A和B两门考试,每门考试得分是不一样的。


# 6个学生
> student <- data.table(id=1:6,name=c('Dan','Mike','Ann','Yang','Li','Kate'));student
   id name
1:  1  Dan
2:  2 Mike
3:  3  Ann
4:  4 Yang
5:  5   Li
6:  6 Kate

# 分别参加A和B两门考试
> score <- data.table(id=1:12,stuId=rep(1:6,2),score=runif(12,60,99),class=c(rep('A',6),rep('B',6)));score
    id stuId    score class
 1:  1     1 89.18497     A
 2:  2     2 61.76987     A
 3:  3     3 74.67598     A
 4:  4     4 64.08165     A
 5:  5     5 85.00035     A
 6:  6     6 95.25072     A
 7:  7     1 81.42813     B
 8:  8     2 82.16083     B
 9:  9     3 69.53405     B
10: 10     4 89.01985     B
11: 11     5 96.77196     B
12: 12     6 97.02833     B

通过学生ID,把学生和考试成绩2个数据集进行连接。


# 设置score数据集,key为stuId
> setkey(score,"stuId")

# 设置student数据集,key为id
> setkey(student,"id")

# 合并两个数据集的数据
> student[score,nomatch=NA,mult="all"]
    id name i.id    score class
 1:  1  Dan    1 89.18497     A
 2:  1  Dan    7 81.42813     B
 3:  2 Mike    2 61.76987     A
 4:  2 Mike    8 82.16083     B
 5:  3  Ann    3 74.67598     A
 6:  3  Ann    9 69.53405     B
 7:  4 Yang    4 64.08165     A
 8:  4 Yang   10 89.01985     B
 9:  5   Li    5 85.00035     A
10:  5   Li   11 96.77196     B
11:  6 Kate    6 95.25072     A
12:  6 Kate   12 97.02833     B

最后我们会看到,两个数据集的结果合并在了一个结果数据集中。这样就完成了,数据连接的操作。从代码的角度来看,1行代码要比用data.frame去拼接方便的多。

3. data.table包性能对比

现在很多时候我们需要处理的数据量是很大的,动辄上百万行甚至上千万行。如果我们要使用R对其进行分析或处理,在不增加硬件的条件下,就需要用一些高性能的数据包进行数据的操作。这里就会发现data.table是非常不错的一个选择。

3.1 data.table和data.frame索引查询性能对比

我们先生成一个稍大数据集,包括2列x和y分别用英文字母进行赋值,100,000,004行,占内存大小1.6G。分别比较data.frame操作和data.table操作的索引查询性能耗时。

使用data.frame创建数据集。


# 清空环境变量
> rm(list=ls())

# 设置大小
> size = ceiling(1e8/26^2)
[1] 147929

# 计算data.frame对象生成的时间 
> t0=system.time(
+   df <- data.frame(x=rep(LETTERS,each=26*size),y=rep(letters,each=size))
+ )

# 打印时间
> t0
用户 系统 流逝 
3.63 0.18 3.80 

# df对象的行数
> nrow(df)
[1] 100000004

# 占用内存
> object.size(df)
1600003336 bytes

# 进行条件查询
> t1=system.time(
+   val1 <- dt[dt$x=="R" & dt$y=="h",]
+ )

# 查询时间
> t1
用户 系统 流逝 
8.53 0.84 9.42 

再使用data.table创建数据集。


# 清空环境变量
> rm(list=ls())

# 设置大小
> size = ceiling(1e8/26^2)
[1] 147929

# 计算data.table对象生成的时间 
> t3=system.time(
+   dt <- data.table(x=rep(LETTERS,each=26*size),y=rep(letters,each=size))
+ )

# 生成对象的时间
> t3
用户 系统 流逝 
3.22 0.39 3.63 

# 对象行数
> nrow(dt)
[1] 100000004

# 占用内存
> object.size(dt)
2000004040 bytes

# 进行条件查询
> t3=system.time(
+ val2 <- dt[x=="R" & y=="h",]
+ )

# 查询时间
> t3
用户 系统 流逝 
6.52 0.26 6.80 

从上面的测试来看,创建对象时,data.table比data.frame显著的高效,而查询效果则并不明显。我们对data.table数据集设置索引,试试有索引查询的效果。


# 设置key索引列为x,y
> setkey(dt,x,y)

# 条件查询
> t4=system.time(
+   val3  <- dt[list("R","h")]
+ )

# 查看时间
> t4
用户 系统 流逝 
0.00 0.00 0.06 

设置索引列后,按索引进行查询,无CPU耗时。震惊了!!

3.2 data.table和data.frame的赋值性能对比

对于赋值操作来说,通常会分为2个动作,先查询再值替换,对于data.frame和data.table都是会按照这个过程来实现的。从上一小节中,可以看到通过索引查询时data.table比data.frame明显的速度要快,对于赋值的操作测试,我们就要最好避免复杂的查询。

对x列值为R的行,对应的y的值进行赋值。首先测试data.frame的计算时间。


> size = 1000000
> df <- data.frame(x=rep(LETTERS,each=size),y=rnorm(26*size))
> system.time(
+   df$y[which(df$x=='R')]<-10
+ )
用户 系统 流逝 
0.75 0.01 0.77 

计算data.table的赋值时间。


> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> system.time(
+   dt[x=='R', y:=10]
+ )
用户 系统 流逝 
0.11 0.00 0.11 
> setkey(dt,x)
> system.time(
+   dt['R', y:=10]
+ )
用户 系统 流逝 
0.01 0.00 0.02 

通过对比data.table和data.frame的赋值测试,有索引的data.table性能优势是非常明显的。我们增大数据量,再做一次赋值测试。


> size = 1000000*5
> df <- data.frame(x=rep(LETTERS,each=size),y=rnorm(26*size))
> system.time(
+   df$y[which(df$x=='R')]<-10
+ )
用户 系统 流逝 
3.22 0.25 3.47 

> rm(list=ls())
> size = 1000000*5
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> setkey(dt,x)
> system.time(
+   dt['R', y:=10]
+ )
用户 系统 流逝 
0.08 0.01 0.08 

对于增加数据量后data.table,要比data.frame的赋值快更多倍。

3.3 data.table和tapply分组计算性能对比

再对比一下data.table处理数据和tapply的分组计算的性能。测试同样地只做一个简单的计算设定,比如,对一个数据集按x列分组对y列求和。


# 设置数据集大小
> size = 100000
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))

# 设置key为x列
> setkey(dt,x)

# 计算按x列分组,对y列的求和时间
> system.time(
+ r1<-dt[,sum(y),by=x]
+ )
用户 系统 流逝 
0.03 0.00 0.03 

# 用tapply实现,计算求和时间
> system.time(
+ r2<-tapply(dt$y,dt$x,sum)
+ )
用户 系统 流逝 
0.25 0.05 0.30 

# 查看数据集大小, 40mb
> object.size(dt)
41602688 bytes

对于40mb左右的数据来说,tapply比data.table要快,那么我增加数据集的大小,给size*10再测试一下。


> size = 100000*10
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> setkey(dt,x)
> val3<-dt[list("R")]
 
> system.time(
+   r1<-dt[,sum(y),by=x]
+ )
用户 系统 流逝 
0.25 0.03 0.28 

> system.time(
+   r2<-tapply(dt$y,dt$x,sum)
+ )
用户 系统 流逝 
2.56 0.36 2.92 

# 400mb数据 
> object.size(dt)
416002688 bytes

对于400mb的数据来说,data.table的计算性能已经明显优于tapply了,再把数据时增加让size*5。


> size = 100000*10*5
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> setkey(dt,x)
 
> system.time(
+     r1<-dt[,sum(y),by=x]
+ )
用户 系统 流逝 
1.50 0.11 1.61 

> system.time(
+     r2<-tapply(dt$y,dt$x,sum)
+ )
 用户  系统  流逝 
13.30  3.58 16.90 
 
# 2G数据
> object.size(dt)
2080002688 bytes

对于2G左右的数据来说,tapply总耗时到了16秒,而data.table为1.6秒,从2个的测试来说,大于400mb数据时CPU耗时是线性的。

把上几组测试数据放到一起,下图所示。

data-table

通过上面的对比,我们发现data.table包比tapply快10倍,比data.frame赋值操作快30倍,比data.frame的索引查询快100倍,绝对是值得花精力去学习的一个包。

赶紧用data.table包去优化你的程序吧!

转载请注明出处:
http://blog.fens.me/r-data-table/

打赏作者