• 粉丝日志首页

用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

打赏作者

EpiModel快速上手传染病模型

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

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

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

前言

在上一篇文章 用R语言解读传染病模型 中,我们已经了解传染病动力学模型的背景知识,并且用R语言手动把基本的传染病模型公式进行了编程实现。对于传染病学科中,已经有专业的团队开发出了专业的模型,所以其实我们不需要手动从头写代码,站在巨人的肩膀前进会让我们进步的更快。

本次将介绍的 EpiModel 包就是一个专业做传染病动力学模型的工具,不仅支持基本的SI, SIS, SIR模型,还支持扩展,自定义开发等酷炫功能,让我们进行传染病的模型开发和模拟,有非常大的帮助。

用R进行跨界建模,打破自己的知识壁垒。世界真奇妙,让我们用R来探索吧。

目录

  1. EpiModel包介绍
  2. 确定性间隔模型的核心函数
  3. SI模型实现
  4. SIS模型实现
  5. SIR模型实现
  6. SEIR模型实现
  7. 基于SI模型实现组间交叉传染
  8. 网络模型实现

1. EpiModel包介绍

EpiModel 提供了用于模拟和分析传染病动力学数学模型的工具,支持的流行病模型类包括确定性隔间模型、随机个体接触模型和随机网络模型。疾病类型包括有人口统计的 SI、SIR 和 SIS 流行病,可实现自定义可扩展的和模拟任意复杂性的流行病模型。EpiModel中 特有的网络模型,基于在 R 的 Statnet 软件套件中实现的时间指数随机图模型 (ERGM) 的统计框架。

官方网站: http://www.epimodel.org/index.html

EpiModel包整体架构

EpiModel包,支持三类流行病模型提供功能:

  • 确定性隔间模型(Deterministic Compartmental Models):模型是确定性的,时间是连续的,解是输入参数和初始条件的固定数学函数,在疾病和人口过渡过程中没有随机变异性。
  • 随机个体接触模型(Stochastic Individual Contact Models):参数是随机抽取,其中控制状态间转移的所有比率和风险都是随机抽取的,它们来自由这些比率或概率总结的分布,包括正态分布、泊松分布和二项分布。时间是离散的,ICMs是在离散时间内模拟的,而DCMs是在连续时间内模拟的。
  • 随机网络模型(Stochastic Network Models): 两个节点之间连边与否不再是确定的事情,而是根据一个概率决定。网络动力学和疾病传播动力学,这两个动态过程是被视为独立的或者相互依赖的。独立的网络模型假设疾病模拟对时间网络的结构没有影响,尽管时间网络的结构肯定会影响疾病。依赖网络模型允许流行病学和人口统计过程影响网络结构。

EpiModel包,默认支持三种传染病类型,同时允许使用者自己扩展模型内核计算,本文介绍的SEIR模型,就是基于自定模块按公式自己开发的。

  • SI:Susceptible-Infective,一种两种状态的疾病,其中存在终生感染而无法恢复。 HIV/AIDS 就是一个例子,尽管在这种情况下,通常将感染阶段建模为单独的隔间。
  • SIR:Susceptible-Infective-Removed,一种三期疾病,感染后可终生恢复并产生免疫力。麻疹就是一个例子,但该疾病的现代模型也需要考虑人群中的疫苗接种模式。
  • SIS:Susceptible-Infective- Susceptible,一种两阶段疾病,其中一个人可能在一生中从易感状态到受感染状态来回转换。例子包括细菌性性传播疾病,如淋病。

EpiModel包框架下,不仅提供了丰富的功能,有很好的扩展性,还支持自定模型的开发,允许多种状态的设计,支持确定计算的和随机的计算,可模拟真实的传染病传播过程。

2. 确定性间隔模型的核心函数

本文中SI,SIS,SIR,SEIR都是基于确定性间隔模型实现的,因此需要重点了解一下确定性间隔模型的核心函数。

主要有4个:

  • param.dcm(),用于设置确定性隔间模型传播参数
  • init.dcm(),用于设置确定性隔间模型的人数相关的初始条件
  • control.dcm(),用于设置确定性隔间模型的控制参数
  • dcm(),加载参数变量,运行模型

查看 param.dcm() 函数定义

> param.dcm
function (inf.prob, inter.eff, inter.start, act.rate, rec.rate,
a.rate, ds.rate, di.rate, dr.rate, inf.prob.g2, act.rate.g2,
rec.rate.g2, a.rate.g2, ds.rate.g2, di.rate.g2, dr.rate.g2,
balance, ...)

参数列表:

查看 init.dcm() 函数定义

> init.dcm
function (s.num, i.num, r.num, s.num.g2, i.num.g2, r.num.g2, ...) 

3. SI模型实现

本文将直接介绍使用EpiModel包的具体实现,关于SI,SIS,SIR模型的具体数据公式和解释,请参考文章:用R语言解读传染病模型

安装和加载EpiModel包

# 安装EpiModel包
> install.packages("EpiModel")

# 加载EpiModel包,并设定运行目录
> library(EpiModel)
> setwd("C:/work/R/covid19")

构建SI模型:初始:日感染率 inf.prob=0.5,日传播率 act.rate=0.5。初始易感染者s.num=500,初始已感染者i.num=1,观察周期nsteps=100。

# 设置模型传播参数
> param <- param.dcm(inf.prob = 0.5, act.rate = 0.5)

# 设置模型人群参数
> init <- init.dcm(s.num = 500, i.num = 1)

# 设置使用SI模型和观察周期
> control <- control.dcm(type = "SI", nsteps = 100)

在param.dcm()函数中,把传染率拆分成了2个值inf.prob和 act.rate,从源代码中这两个参数相乘就是日感染率,这里是给我们提供了一个额外的调解参数,并没有改变模型本身的算法。

运行模型,并查看模型计算结果。查看mod变量,会打印SI模型所有运行参数。

# 运行模型
> mod <- dcm(param, init, control)
> mod
EpiModel Simulation
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SI
No. runs: 1
No. time steps: 100
No. groups: 1

Model Parameters
-----------------------
inf.prob = 0.5
act.rate = 0.5

Model Output
-----------------------
Variables: s.num i.num si.flow num

模型结果,可视化输出,直接使用plot()函数即可。

> plot(mod)


上图中,蓝色线为易感染人数,红色线为已感染人数,从图中我们能看到人群的变化过程。

查看模型结果数据集,前10条计算结果,可以直接把mod 转型为data.frame类型,来查看每一个周期的运行结果。

> head(as.data.frame(mod),10)
   time    s.num    i.num   si.flow num
1     1 500.0000 1.000000 0.2832895 501
2     2 499.7167 1.283290 0.3632782 501
3     3 499.3534 1.646568 0.4656817 501
4     4 498.8878 2.112249 0.5966712 501
5     5 498.2911 2.708921 0.7640466 501
6     6 497.5270 3.472967 0.9776201 501
7     7 496.5494 4.450587 1.2496619 501
8     8 495.2998 5.700249 1.5953944 501
9     9 493.7044 7.295644 2.0335064 501
10   10 491.6708 9.329150 2.5866250 501

如果想查看第20个观察周期的传播率、易感染者人数和已感染者人数时,可以使用comp_plot()进行画图,把这一时刻的传播关系可视化展示。

> comp_plot(mod, at = 20, digits = 2)

如果想查看第20个观察周期的具体的计算数值时,可以使用summary()函数,指定at=20来看模型的明细数据。

> summary(mod, at = 20)
EpiModel Summary
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SI
No. runs: 1
No. time steps:
No. groups: 1

Model Statistics
------------------------------
Time: 20	 Run: 1 
------------------------------ 
                n    pct
Suscept.  406.938  0.812
Infect.    94.062  0.188
Total     501.000  1.000
S -> I     20.601     NA
------------------------------ 

通过EpiMdoel包来计算SI模型的模拟运行,是非常方便的,省去了我们原来自定义微分方程的过程,而且EpiModel还集成是可视化图、传播链图,参数多样化等高级特性,让我们做传染病的模拟,仅需要几行代码就可以实现。

4. SIR模型实现

通过 SIR模型,进行人口分析。在易感-感染-恢复 (SIR) 模型中,受感染的个体从疾病转变为终生恢复状态,在这种状态下他们将不再易感。

初始:日感染率inf.prob=0.2,日传播率 act.rate=1,恢复率rec.rate=1/20, 出生率a.rate=1/95,死亡率ds.rate=1/100,疾病死亡率 di.rate=1/80, 恢复率dr.rate=1/100。易感染者s.num=1000,已感染者i.num=1,治愈r.num=0,观察周期nsteps=500。


# 设置模型传播参数
> param <- param.dcm(inf.prob = 0.2, act.rate = 1, rec.rate = 1/20,
+                    a.rate = 1/95, ds.rate = 1/100, di.rate = 1/80,  r.rate = 1/100)

# 设置模型人群参数
> init <- init.dcm(s.num = 1000, i.num = 1, r.num = 0)

# 设置使用SIR模型、观察周期、间隔步长
> control <- control.dcm(type = "SIR", nsteps = 500, dt = 0.5)

运行模型,并查看模型计算结果

# 运行模型
> mod <- dcm(param, init, control)
> mod
EpiModel Simulation
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SIR
No. runs: 1
No. time steps: 500
No. groups: 1

Model Parameters
-----------------------
inf.prob = 0.2
act.rate = 1
rec.rate = 0.05
a.rate = 0.01052632
ds.rate = 0.01
di.rate = 0.0125
dr.rate = 0.01

Model Output
-----------------------
Variables: s.num i.num r.num si.flow ir.flow a.flow 
ds.flow di.flow dr.flow num

模型结果,可视化输出

> par(mar = c(3.2, 3, 2, 1), mgp = c(2, 1, 0), mfrow = c(1, 2))
> plot(mod, popfrac = FALSE, alpha = 0.5, lwd = 4, main = "Compartment Sizes")
> plot(mod, y = "si.flow", lwd = 4, col = "firebrick", main = "Disease Incidence", legend = "n")

左图为SIR模型中三种人群的人群,蓝色线s.num易感染人数,红色约i.num为已感染人数,绿色r.num为已治愈人数,由于有了一定的治愈率dr.rate,因此在100周期后3条线基本稳定。右图为发病率,从易感人群到已感染的传播过程。

如果我们想查看在第49个观察周期的传播率数值,易感染者人数,已感染者人数,恢复者人数,以及各状态之间的传播数值。

> par(mfrow = c(1, 1))
> comp_plot(mod, at = 49, digits = 1)

查看前10条结果,和后10条结果。观察num数量,在SIR模型中,总人数增加了,因为出生率值高。

# 查看前10行
> head(as.data.frame(mod))
  time    s.num    i.num      r.num   si.flow    ir.flow   a.flow  ds.flow     di.flow       dr.flow      num
1  1.0 1000.000 1.000000 0.00000000 0.1034038 0.02587806 5.269111 5.000416 0.006469515 0.00006384888 1001.000
2  1.5 1000.165 1.071056 0.02581421 0.1107397 0.02771672 5.270491 5.001226 0.006929180 0.00019713443 1001.262
3  2.0 1000.324 1.147150 0.05333379 0.1185940 0.02968571 5.271870 5.002000 0.007421427 0.00033924715 1001.524
4  2.5 1000.475 1.228637 0.08268026 0.1270029 0.03179423 5.273250 5.002738 0.007948557 0.00049081575 1001.786
5  3.0 1000.619 1.315897 0.11398367 0.1360055 0.03405211 5.274629 5.003435 0.008513027 0.00065251323 1002.048
6  3.5 1000.754 1.409337 0.14738326 0.1456431 0.03646987 5.276008 5.004088 0.009117467 0.00082505998 1002.311

# 查看后10行
> tail(as.data.frame(mod)) 
     time    s.num    i.num    r.num  si.flow  ir.flow   a.flow  ds.flow   di.flow  dr.flow      num
994 497.5 348.0338 129.4081 632.5081 4.057928 3.235373 5.842185 1.740279 0.8088432 3.162722 1109.950
995 498.0 348.0778 129.4218 632.5807 4.058395 3.235716 5.842871 1.740499 0.8089290 3.163085 1110.080
996 498.5 348.1218 129.4355 632.6533 4.058862 3.236060 5.843557 1.740719 0.8090151 3.163448 1110.211
997 499.0 348.1658 129.4493 632.7260 4.059331 3.236405 5.844243 1.740939 0.8091014 3.163811 1110.341
998 499.5 348.2097 129.4631 632.7985 4.059801 3.236752 5.844930 1.741159 0.8091879 3.164174 1110.471
999 500.0 348.2537 129.4770 632.8711 4.059801 3.236752 5.844930 1.741159 0.8091879 3.164174 1110.602

4. SIS模型实现

用SIS模型进行敏感性分析,在易感-感染-恢复 (SIS) 模型中,我们在一个封闭的群体中对此进行,感染力和恢复等式相互反映,因为个体在状态之间来回流动。

建模初始:日感染率inf.prob=(0.1 0.1 0.2 0.1 0.2 0.2),日传播率 act.rate=(0.2 0.4 0.4 0.6 0.6 1),日治愈率 rec.rate=0.02。初始易感染者s.num=500,初始已感染者i.num=1,观察周期nsteps=350。

> inf<-c(0.1, 0.1, 0.2, 0.1, 0.2,0.2)
> act<-c(0.2, 0.4, 0.4, 0.6, 0.6,1)

# 设置模型传播参数
> param <- param.dcm(inf.prob = inf, act.rate = act, rec.rate = 0.02)

# 设置模型人群参数
> init <- init.dcm(s.num = 500, i.num = 1)

# 设置使用SIS模型和观察周期
> control <- control.dcm(type = "SIS", nsteps = 350)

我们同时给出了6组数据,把6组数据的参数进行微调,可一次性得出结果。运行模型,并查看模型计算结果

# 运行模型
> mod <- dcm(param, init, control)
> mod
EpiModel Simulation
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SIS
No. runs: 6
No. time steps: 350
No. groups: 1

Model Parameters
-----------------------
inf.prob = 0.1 0.1 0.2 0.1 0.2 0.2
act.rate = 0.2 0.4 0.4 0.6 0.6 1
rec.rate = 0.02

Model Output
-----------------------
Variables: s.num i.num si.flow is.flow num

模型结果,可视化输出。同时输出6组结果,并画在一张图中,进行对比。

> par(mfrow = c(1,2), mar = c(3.2,3,2.5,1))
> plot(mod, alpha = 1, main = "Disease Prevalence",legend = "full") #已感染者数量
> plot(mod, y = "si.flow", col = "Greens", alpha = 0.8, 
+      main = "Disease Incidence",legend = "full")

左图中,run1,run2,run3,run4,run5,run6分别表示6次运行结果,当inf.prob * act.prob 越大时,疾病从易感染者向易感染者传播越快、曲线越陡,因此run6是最高的曲线。从图中run1几乎没有传播率,是因为治愈率与传播率相等,相互低效了。右图中,为6次模拟的发病率。

SIS在计算时,根据不同的参数输入,我们共跑了6次模型,如果想查看其中某次的计算结果,可以使用run参数。

# 查看第1次运行结果的前6行输出
> head(as.data.frame(mod, run = 1)) #提取输出
  run time    s.num     i.num    si.flow    is.flow num
1   1    1 500.0000 1.0000000 0.01995968 0.01999960 501
2   1    2 500.0000 0.9999601 0.01995889 0.01999880 501
3   1    3 500.0001 0.9999202 0.01995809 0.01999800 501
4   1    4 500.0001 0.9998803 0.01995730 0.01999721 501
5   1    5 500.0002 0.9998403 0.01995650 0.01999641 501
6   1    6 500.0002 0.9998004 0.01995571 0.01999561 501

# 查看第6次运行结果的后6行输出
> tail(as.data.frame(mod, run = 6)) #提取输出
    run time s.num i.num si.flow is.flow num
345   6  345  50.1 450.9   9.018   9.018 501
346   6  346  50.1 450.9   9.018   9.018 501
347   6  347  50.1 450.9   9.018   9.018 501
348   6  348  50.1 450.9   9.018   9.018 501
349   6  349  50.1 450.9   9.018   9.018 501
350   6  350  50.1 450.9   9.018   9.018 501

5. SEIR模型实现

易感-潜伏-感染-恢复 (SEIR) 模型,扩展了 SIR 模型以包括潜伏期但非传染性的类别,模型虑了开放人群中易感人群、潜伏人群、感染者的比例。适用于存在易感者、潜伏者、患病者和康复者四类人群,有潜伏期、治愈后获得终身免疫的疾病,如带状疱疹。

  • S (Susceptible),易感染者,指缺乏免疫能力健康人,与感染者接触后容易受到感染
  • E (Exposed),潜伏者 ,指接触过感染者但不存在传染性的人,可用于存在潜伏期的传染病
  • I (Infectious),已感染者,指有传染性的病人,可以传播给 S,将其变为 E 或 I
  • R (Recovered),康复者,指病愈后具有免疫力的人,如是终身免疫性传染病,则不可被重新变为 S 、E 或 I ,如果免疫期有限,就可以重新变为 S 类,进而被感染。

自定义SEIR模型的内核函数:

> SEIR <- function(t, t0, parms) {
+   with(as.list(c(t0, parms)), {
+     
+     # Population size
+     num <- s.num + e.num + i.num + r.num
+     
+     # Effective contact rate and FOI from a rearrangement of Beta * c * D
+     ce <- R0 / i.dur
+     lambda <- ce * i.num/num
+     
+     dS <- -lambda*s.num
+     dE <- lambda*s.num - (1/e.dur)*e.num
+     dI <- (1/e.dur)*e.num - (1 - cfr)*(1/i.dur)*i.num - cfr*(1/i.dur)*i.num
+     dR <- (1 - cfr)*(1/i.dur)*i.num
+     
+     # Compartments and flows are part of the derivative vector
+     # Other calculations to be output are outside the vector, but within the containing list
+     list(c(dS, dE, dI, dR, 
+            se.flow = lambda * s.num,
+            ei.flow = (1/e.dur) * e.num,
+            ir.flow = (1 - cfr)*(1/i.dur) * i.num,
+            d.flow = cfr*(1/i.dur)*i.num),
+          num = num,
+          i.prev = i.num / num,
+          ei.prev = (e.num + i.num)/num)
+   })
+ }

建模初始:初始繁殖数R0 = 1.9, 潜伏状态的持续时间 e.dur = 10, 染状态的持续时间i.dur = 14, 病死率cfr = c(0.5, 0.7, 0.9)。初始易感染者s.num=1000*1000,初始潜伏者人数e.num=10,初始已感染者i.num=0,初始治愈者r.num=0。观察周期nsteps=350,设置自定模型函数SEIR。

# 设置模型传播参数
> param <- param.dcm(R0 = 1.9, e.dur = 10, i.dur = 14, cfr = c(0.5, 0.7, 0.9))

# 设置模型人群参数
> init <- init.dcm(s.num = 1e6, e.num = 10, i.num = 0, r.num = 0,
+                  se.flow = 0, ei.flow = 0, ir.flow = 0, d.flow = 0)

# 设置使用SEIR模型和观察周期
> control <- control.dcm(nsteps = 500, dt = 1, new.mod = SEIR)

运行模型,并查看模型计算结果


> mod <- dcm(param, init, control)
> mod
EpiModel Simulation
=======================
Model class: dcm

Simulation Summary
-----------------------
No. runs: 3
No. time steps: 500

Model Parameters
-----------------------
R0 = 1.9
e.dur = 10
i.dur = 14
cfr = 0.5 0.7 0.9

Model Output
-----------------------
Variables: s.num e.num i.num r.num se.flow ei.flow ir.flow 
d.flow num i.prev ei.prev

模型结果,可视化输出。通过可视化可分别描述4种状态的人群分布曲线。

par(mfrow = c(2, 2))
plot(mod, y = "i.num", run = 2, main = "患病率")
plot(mod, y = "se.flow", run = 2, main = "发病率")
plot(mod, y = "i.num", main = "感染人数")
plot(mod, y = "i.prev", main = "感染百分比", ylim = c(0, 0.5), legend = "full")

6. 基于SI模型实现组间交叉传染

EpiMoel除了能做单模型的计算,还能进行分组,做组间模型的传播实验。

我们接下来,基于SI模型实现组间交叉传播的模拟,在当前结构的两组模型中,组之间的混合纯粹是异质的:一个组只与另一组有联系,没有组内联系。

比如,在异性恋疾病传播,其中群体代表性别,我们可以将第 1 组建模为女性,将第 2 组建模为男性。

建模初始:日感染率女性inf.prob=0.4,男性inf.prob.g2=0.1。易感染者女性s.num=500男性s.num.g2=500,已感染者女性i.num=1,男性i.num.g2=0 ,观察周期nsteps=500。

# 设置模型传播参数
> param <- param.dcm(inf.prob = 0.4,  inf.prob.g2 = 0.1, act.rate = 0.25, 
+ balance = "g1",a.rate = 1/100, a.rate.g2 = NA, ds.rate = 1/100, 
+ ds.rate.g2 = 1/100, di.rate = 1/50, di.rate.g2 = 1/50)

# 设置模型人群参数
> init <- init.dcm(s.num = 500, i.num = 1, s.num.g2 = 500, i.num.g2 = 0)

# 设置使用SI模型和观察周期
> control <- control.dcm(type = "SI", nsteps = 500)

运行模型,并查看模型计算结果


# 运行模型
> mod <- dcm(param, init, control)
> mod
EpiModel Simulation
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SI
No. runs: 1
No. time steps: 500
No. groups: 2

Model Parameters
-----------------------
inf.prob = 0.4
act.rate = 0.25
a.rate = 0.01
ds.rate = 0.01
di.rate = 0.02
inf.prob.g2 = 0.1
a.rate.g2 = NA
ds.rate.g2 = 0.01
di.rate.g2 = 0.02
balance = g1

Model Output
-----------------------
Variables: s.num i.num s.num.g2 i.num.g2 si.flow a.flow 
ds.flow di.flow si.flow.g2 a.flow.g2 ds.flow.g2 di.flow.g2 
num num.g2

模型结果,可视化输出

> par(mfrow = c(1, 1))
> plot(mod)

上图中,实线为女性组,虚线为男性组。

查看第20个观察周期的数据值

> summary(mod, at = 20)
EpiModel Summary
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SI
No. runs: 1
No. time steps:
No. groups: 2

Model Statistics
------------------------------------------------------------
Time: 20	 Run: 1 
------------------------------------------------------------ 
                   n:g1  pct:g1     n:g2  pct:g2
Suscept.        499.803   0.998  499.746   0.999
Infect.           1.016   0.002    0.376   0.001
Total           500.819   1.000  500.122   1.000
S -> I            0.038      NA    0.026      NA
Arrival ->        5.008      NA    5.008      NA
S Departure ->    4.998      NA    4.997      NA
I Departure ->    0.021      NA    0.008      NA
------------------------------------------------------------ 

5. 网格模型实现

使用 EpiModel 的基本网络模型,是基于时间指数随机图模型的框架构建任意复杂的联系或伙伴关系结构,这些结构会随着时间的推移而形成和消散。

官方例子:血清监测:疾病状况影响伴侣选择——人口统计过程(出生、死亡和迁移)——接触网络过程必须适应人口规模和组成的变化。对于独立模型,在每次疫情模拟开始时对全动态网络进行模拟。对于相关模型,在每个时间步长的情况下,网络被重新模拟成变化的种群大小和变化的节点属性的函数。

设置网络模型

> nw <-network.initialize(n = 1000, directed = FALSE)
> nw <-set.vertex.attribute(nw, "race", rep(0:1, each = 500))
> formation <-~edges + nodefactor("race") + nodematch("race") +concurrent
> target.stats <- c(250, 375, 225, 100)
> coef.diss <-dissolution_coefs(dissolution = ~offset(edges), duration = 25)
> est1 <- netest(nw, formation,target.stats, coef.diss, edapprox = TRUE)
Starting maximum pseudolikelihood estimation (MPLE):
Evaluating the predictor and response matrix.
Maximizing the pseudolikelihood.
Finished MPLE.
Starting Monte Carlo maximum likelihood estimation (MCMLE):
Iteration 1 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.1944.
Estimating equations are not within tolerance region.
Iteration 2 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.0656.
Convergence test p-value: 0.4882. Not converged with 99% confidence; increasing sample size.
Iteration 3 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.1167.
Estimating equations are not within tolerance region.
Iteration 4 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.0288.
Convergence test p-value: 0.0892. Not converged with 99% confidence; increasing sample size.
Iteration 5 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.0142.
Convergence test p-value: 0.0555. Not converged with 99% confidence; increasing sample size.
Iteration 6 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.0154.
Convergence test p-value: 0.0138. Not converged with 99% confidence; increasing sample size.
Iteration 7 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.0047.
Convergence test p-value: 0.0002. Converged with 99% confidence.
Finished MCMLE.
This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the
mcmc.diagnostics() function.

进行模型诊断

> #模型诊断
> dx <- netdx(est1, nsims = 5, nsteps =500,
+             nwstats.formula = ~edges + nodefactor("race", base = 0) + nodematch("race") + concurrent)

Network Diagnostics
-----------------------
- Simulating 5 networks
  |In term ‘nodefactor’ in package ‘ergm’: Argument ‘base’ has been superseded by ‘levels’, and it is recommended to use the latter.  Note that its interpretation may be different.
*****|
- Calculating formation statistics
- Calculating duration statistics
- Calculating dissolution statistics

设置模型参数

> param <- param.net(inf.prob = 0.1,act.rate = 5, rec.rate = 0.02)
> status.vector <- c(rbinom(500, 1, 0.1),rep(0, 500)) # 二项分布
> status.vector <- ifelse(status.vector ==1, "i", "s")
> init <- init.net(status.vector =status.vector)
> control <- control.net(type ="SIS", nsteps = 500, nsims = 10, epi.by = "race")

模型构建,进行各个阶段在10个周期的传播。(时间较长)

> sim1 <- netsim(est1, param, init,control)
Epidemic Simulation
----------------------------
Simulation: 10/10
Timestep: 500/500
Prevalence: 459
Population Size: 1000
----------------------------

# 查看模型运行状态和参数
> sim1
EpiModel Simulation
=======================
Model class: netsim

Simulation Summary
-----------------------
Model type: SIS
No. simulations: 10
No. time steps: 500
No. NW groups: 1

Fixed Parameters
---------------------------
inf.prob = 0.1
act.rate = 5
rec.rate = 0.02
groups = 1

Model Output
-----------------------
Variables: s.num s.num.race0 s.num.race1 i.num i.num.race0 
i.num.race1 num num.race0 num.race1 si.flow is.flow
Networks: sim1 ... sim10
Transmissions: sim1 ... sim10

Formation Diagnostics
----------------------- 
                  Target Sim Mean Pct Diff Sim SD
edges                250  250.084    0.034 14.896
nodefactor.race.1    375  375.007    0.002 25.379
nodematch.race       225  225.192    0.085 14.044
concurrent           100   99.629   -0.371 10.772


Dissolution Diagnostics
----------------------- 
               Target Sim Mean Pct Diff Sim SD
Edge Duration   25.00   24.743   -1.026  0.295
Pct Edges Diss   0.04    0.040    0.614  0.001

查看模型运行结果


> head(as.data.frame(sim1), 10)
   sim time s.num s.num.race0 s.num.race1 i.num i.num.race0 i.num.race1  num num.race0 num.race1 si.flow is.flow
1    1    1   938         438         500    62          62           0 1000       500       500      NA      NA
2    1    2   933         434         499    67          66           1 1000       500       500       6       1
3    1    3   933         435         498    67          65           2 1000       500       500       2       2
4    1    4   933         438         495    67          62           5 1000       500       500       4       4
5    1    5   927         436         491    73          64           9 1000       500       500       7       1
6    1    6   928         438         490    72          62          10 1000       500       500       3       4
7    1    7   926         437         489    74          63          11 1000       500       500       3       1
8    1    8   924         438         486    76          62          14 1000       500       500       3       1
9    1    9   925         440         485    75          60          15 1000       500       500       2       3
10   1   10   926         440         486    74          60          14 1000       500       500       1       2

获得网络状态

# 获得网络状态
> nw <- get_network(sim1, sim = 1)

#获得每个时间点的时间状态
get_transmat(sim1, sim = 1)

画图输出2个时间点的网络状态:

> par(mfrow = c(1,2), mar = c(0,0,1,0))
> plot(sim1, type = "network", at =1, col.status = TRUE,
+      main = "Prevalence at t1")
> plot(sim1, type = "network", at =500, col.status = TRUE,
+      main = "Prevalence at t500")

上图中左边为第1个周期的状态,右边为第500个周期的状态。

查看第500个周期的数值


> summary(sim1,at=500)

EpiModel Summary
=======================
Model class: netsim

Simulation Details
-----------------------
Model type: SIS
No. simulations: 10
No. time steps: 500
No. NW groups: 1

Model Statistics
------------------------------
Time: 500 
------------------------------ 
            mean      sd    pct
Suscept.   529.4  13.066  0.529
Infect.    470.6  13.066  0.471
Total     1000.0   0.000  1.000
S -> I      10.8   2.821     NA
I -> S       8.5   2.915     NA
------------------------------ 

通过使用EpiMdoel包,让我们做传染病模型的传播模拟,快速、高效。在掌握了,传染病的基础知识后,少量的代码就可以让我们做出专业性很强的结果,确实是站在巨人的肩膀,让知识传播更快、更高效。

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

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

打赏作者

用R语言解读传染病模型

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

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

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

前言

传染病学是人们在不断地同危害人类健康严重的疾病作斗争中发展起来的,新冠疫情已经持续三年了,近期在北京和上海两个超大城市,还在不停的变种和传播,居家办公,每天核酸,都已经变成了生活的常态。

好奇心驱使我要研究一下传染病,为什么疫情会有这么大的传染性,如果不做防控会怎么样。结合传染病学的专业知识,通过R语言的技术手段,彻底明白了传染病传播的一些底层原理。

由于本人非传染病专业,文中关于传染病的专业内容大都从网上汇集和再整理,当出现本文描述与教课书不符时,请以专业教材为准。

目录

  1. 传染病动力学介绍
  2. 自由增长模型
  3. SI模型
  4. SIS模型
  5. SIR模型

1. 传染病动力学介绍

在全球抗击新冠肺炎疫情的过程中,传染病动力学模型发挥了巨大的作用。传染病动模型,根据传染病的传播速度不同,空间范围各异,传播途径多样,动力学机理等各种因素,对传染病模型按照传染病的类型划分为 SI,SIR,SIRS,SEIR 模型。

对传染病的研究方法主要有四种:描述性研究、分析性研究、实验性研究和理论性研究。传染病动力学是对传染病进行理论性定量研究的一种重要方法。它是根据种群生长的特性,疾病的发生及在种群内的传播、发展规律,以及与之有关的社会等因素,建立能反映传染病动力学特性的数学模型。通过对模型动力学性态的定性、定量分析和数值模拟,来显示疾病发展程,寻找对期预防和控制的最优策略。

2. 自由增长模型

自由增长模型,是一种条件最简单的传染病模型,模型整体上呈现指数增长趋势,对传染病传播率ρ非常敏感。模型适合于古代医疗条件不发达、不懂得对病人进行防疫隔离时,发生恶性瘟疫的情形。

模型假设:

  • I(t)为t时刻的感染者人数,I(t)是连续、可微函数。
  • ρ是常数,为每个病人每天有效接触的人数,即足以使人致病的接触, ρ >0

模型构建:

  • 从t到t+Δt的时间里,病人增加的人数为
I(t+Δt) – x(t) = ρ * I(t)* Δt

模型求解:

  • 当t=0时,初始病人人数为 I0 得,

𝑑𝐼/𝑑𝑡= ρ * I
I(0) = I0

R语言代码实现


# 加载程序包
> library(deSolve)  # 用于计算积分
> library(magrittr) 
> library(ggplot2)
> library(reshape2)
> setwd("C:/work/R/covid19")

定义2个公共函数,draw()函数用于画图,rbinddf()函数用于list转型为data.frame。


# 画图函数
> draw<-function(sdf,title="SIR模型",ylab="感染人数",xlab="时间周期",ylog=FALSE){
+   g<-ggplot(data = sdf, mapping = aes(x=time,y=value,colour=type))
+   g<-g+geom_line() + geom_point()                                   #绘制线图和点图
+   g<-g+scale_shape_manual(values = c(21,23))                        #自定义点形状
+   if(ylog) g<-g+scale_y_log10()
+   g<-g+ggtitle(title)+xlab(xlab)+ylab(ylab)
+   g
+ }

# 按行合并
> rbinddf<-function(ldf=list(),names=NULL){
+   n <- length(ldf)
+   res <- NULL
+   if(n>0){
+     for (i in seq(n)) {
+       if(!is.null(names)){
+         type <- names[i]
+       }
+       res <- rbind(res, data.frame(ldf[[i]],type=type))
+     }
+   }
+   return(res)
+ }

定义自由增长模型内核函数


# case1 自由增长模型
> calc1 <- function(rio = 0.3,         # 传染率系数,每个病人每天传染的人数
+                   times = 0:20,      # 病情发展时间
+                   I = 1) {           # 已感染者1人
+   # 公式
+   s_equations <- function(times, vars, params) {
+     with(as.list(c(vars, params)), {
+       dI <- rio * I
+       return(list(dI))
+     })
+   }
+   
+   # 解微分方程
+   ode(
+     func = s_equations,
+     y = c(I = I),
+     times = times,
+     parms = c(rio = rio)
+   ) %>%  as.data.frame
+ }

分别计算当ρ(rio) = 0.3, 0.25, 0.2 时的感染者人员。


> s1 <- calc1(rio = 0.3)
> s2 <- calc1(rio = 0.25)
> s3 <- calc1(rio = 0.2)

# 合并3组结果
> sdf1<-rbinddf(list(s1,s2,s3),names=c("s1","s2","s3"))
> head(sdf1)         # 查看数据集
  time        I type
1    0 1.000000   s1
2    1 1.349861   s1
3    2 1.822121   s1
4    3 2.459609   s1
5    4 3.320123   s1
6    5 4.481699   s1

把模拟进行可视化展示。


> names(sdf1)<-c("time","value","type")
> draw(sdf1,title="自由增长模型",ylab="感染人数",xlab="时间周期")

输出如图所示

自由增长模型,并不符合实际情况。在病人有效接触的人群中,有健康人也有病人,而其中只有健康人才可以被传染为病人,所以我们接下来介绍的SI模型,就区别2类人群。

3. SI模型

SI模型,适用于只有易感者和患病者两类人群,且不会反复发作的疾病。人群分为易感染者(Susceptible)和已感染者(Infective)两类(取两个词的首字母,称之为SI模型)。

易感者与患病者有效接触即被感染,变为患病者,无潜伏期、无治愈情况、无免疫力。SI模型对传染病传播率ρ非常敏感,而且只要我们把每个病人、每天接触的人数有效降低,传染病的传染速度就会变得非常慢。只要防疫力度够大,控制住传染病是完全可能的。这种模型比较适合,被传染后无法恢复健康的,比如HIV等情形。

模型假设:

  • I(t)为t时刻的感染者人数,I(t)是连续、可微函数。
  • ρ是常数,为每个病人每天有效接触的人数,即足以使人致病的接触, ρ >0
  • 在疾病传播期内所考察地区的总人数不变,不考虑生死和迁移。
  • 当病人与健康者有效接触时,使健康者受感染变为病人

模型构建:
已感染人数为 N * i(t),所有已感染人每天有效接触人数为 ρ * N * i(t),其中被传染的易感染人的比例 s(t),因此每天增加的感染者人数为 ρ * N * i(t) * s(t)。

  • N:总人数,不变
  • s(t):t时刻易感染者在总人数中所占比例。
  • i(t):t时刻已感染者在总人数中所占比例。
  • ρ是常数,每个病人每天有效接触的平均人数是常数,称为日接触率, ρ >0

从t到t+Δt的时间里,病人增加的人数为

N * (i(t+Δt)– i(t)) = ρ *  i(t) * s(t) * Δt

模型求解:

  • N总人数不变,因此 s(t) + i(t) = 1
  • 当t=0时,初始病人占比为 i0 得,

𝑑𝑖/𝑑𝑡= ρ * i * (1-i)
i(0) = 𝑖0

编写模型内核函数。


> calc2 <- function(rio = 0.3,                # 接触率,传染率系数
+                  times = 0:100,             # 病情发展时间
+                  i = 0.000001) {            # 已感染者,初始占比:百万分之一
+   # 公式
+   si_equations <- function(time, vars, params) {
+     with(as.list(c(vars, params)), {
+       di <- rio * i * (1 - i)
+       return(list(di))
+     })
+   }
+   
+   # 解微分方程
+   ode(
+     func = si_equations,  
+     y = c(i = i),         
+     times = times,        
+     parms = c(rio = rio)
+   ) %>%  as.data.frame
+ }

分别计算当传播率ρ(rio)为0.3,0.25,0.2,0.15,0.12时的结果,已感染者比例𝑖0 =0.000001(百万分之一)


> si1 <- calc2(rio = 0.3)
> si2 <- calc2(rio = 0.25)
> si3 <- calc2(rio = 0.2)
> si4 <- calc2(rio = 0.15)
> si5 <- calc2(rio = 0.12)
> sdf2<-rbinddf(list(si1,si2,si3,si4,si5),
+              names=c("si1","si2","si3","si4","si5"))
> names(sdf2)<-c("time","value","type")

# 画图
> draw(sdf2,title="SI模型",ylab="感染人数占比",xlab="时间周期")

4. SIS模型

人群分为易感染者(Susceptible)和已感染者(Infective)两类,已感染者(Infective)可以被治愈,变成易感染者(Susceptible)(称之为SIS模型)。SIS模型与SI模型的差异,在于SIS模型假设已感染者(Infective)可以被治愈,重新变成易感染者(Susceptible),比如季节性流感等。

1、传染的速度,取决于病人的传染速度,与病人的治愈速度之间的相对水平,如果病人的治愈速度,大于病人的传染速度,即μ>ρ,那么该传染病最终会消失;
2、即便病人的传染速度高于治愈速度,最终也只有一部分人群会被感染;
3、最终感染的人群比例,为1-1/sigma,sigma = ρ / μ ,sigma是整个传染期内每个病人的有效接触人数,可以理解为病人在整个生病期内,接触的总人数。

模型假设

  • 在疾病传播期内所考察地区的总人数不变,不考虑生死和迁移。
  • 当病人与健康者有效接触时,使健康者受感染变为病人。
  • 人被全部治愈,所需要的天数为1/μ,为传染病的平均传染期。

模型构建:
已感染人数为 N * i(t),所有已感染人每天有效接触人数为 ρ * N * i(t),其中被传染的易感染人的比例 s(t),因此每天增加的感染者人数为 ρ * N * i(t) * s(t)。被治愈的人数为 μ * N * i(t)。

  • N:总人数,不变
  • s(t):t时刻易感染者在总人数中所占比例。
  • i(t):t时刻已感染者在总人数中所占比例。
  • ρ是常数,每个病人每天有效接触的平均人数是常数,为日接触率, ρ >0
  • μ是常数,每天被治愈的病人数占病人总数的比例,为日治愈率

从t到t+Δt的时间里,病人增加的人数为

N * (i(t+Δt)– i(t)) = ρ *  i(t) * s(t) * Δt – μ * N * i(t) * Δt

模型求解:

  • N总人数不变,因此 s(t) + i(t) = 1
  • 当t=0时,初始病人占比为 i0 得,

𝑑𝑖/𝑑𝑡= ρ * i * (1-i) – μ * i
i(0) = 𝑖0

编写模型内核函数。


> calc3 <- function(rio = 0.3,                 # 接触率,传染率系数
+                   mu=0.1,                    # 治愈率
+                   times = 0:150,             # 病情发展时间
+                   i = 0.000001) {            # 已感染者初始占比:百万分之一
+   sis_equations<-function(time,vars,params){
+     with(as.list(c(vars,params)),{
+       di<-rio*i*(1-i)-mu*i
+       return(list(c(di)))
+     })
+   }
+   
+   ode(
+     func=sis_equations,
+     y=c(i=i),
+     times=times,
+     parms=c(rio=rio,mu=mu)
+   ) %>% as.data.frame()
+ }

分别日接触率 ρ(rio) = 0.3、0.25、0.2、0.15、0.12,日治愈率 μ =0.1,已感染者比例 𝑖0 =0.000001(百万分之一)


> sis1 <- calc3(rio = 0.27,mu=0.1)
> sis2 <- calc3(rio = 0.25,mu=0.1)
> sis3 <- calc3(rio = 0.20,mu=0.1)
> sis4 <- calc3(rio = 0.15,mu=0.1)
> sis5 <- calc3(rio = 0.12,mu=0.1)
> sdf<-rbinddf(list(sis1,sis2,sis3,sis4,sis5),
+              names=c("sis1","sis2","sis3","sis4","sis5"))
> names(sdf)<-c("time","value","type")
> draw(sdf,title="SIS模型",ylab="感染人数占比",xlab="时间周期")

5. SIR模型

人群分为易感染者(Susceptible)和已感染者(Infective)和恢复者(Removed)三类(称之为SIR模型)。SIR模型与SIS模型的差异,在于SIR模型假设已感染者(Infective)被治愈后会具备免疫力,不会被感染,从而成为恢复者(Removed),SIR模型已经初步接近实际的传染病模型。

模型假设:

  • 在疾病传播期内所考察地区的总人数不变,不考虑生死和迁移。
  • 当病人与健康者有效接触时,使健康者受感染变为病人。
  • 人被全部治愈,所需要的天数为1/μ,为传染病的平均传染期。

模型构建:
已感染人数为 N * i(t),已感染人每天有效接触人数为 ρ * I(t),被传染的易感染人的比例 S(t)/N,因此每天增加的感染者人数为 ρ * I(t) * S(t)/N,被治愈的人数为 μ * I(t)。

  • N:总人数,不变
  • S(t):t时刻易感染者在总人数中所占比例。
  • I(t):t时刻已感染者在总人数中所占比例。
  • R(t): t时刻恢复者在总人数中所占比例。
  • ρ是常数,每个病人每天有效接触的平均人数是常数,为日接触率, ρ >0
  • μ是常数,每天被治愈的病人数占病人总数的比例,为日治愈率

从t到t+Δt的时间里,易感染者人数,已感染者人数,恢复者人数计算:


# 易感染者人数计算
S(t+ Δt) – S(t) = - ρ * I(t) * Δt/N   

# 已感染者人数计算
I(t+ Δt) – I(t) = ρ * I(t) * S(t) * Δt/N - μ * I(t) * Δt

# 恢复者人数计算
R(t+ Δt) – R(t) = μ * I(t) * Δt

模型求解:

  • N总人数不变,因此 s(t) + i(t) = 1
  • 当t=0时,初始病人占比为 i0 得,

𝑑𝑆/𝑑𝑡= -ρ * I * S/N
𝑑𝐼/𝑑𝑡= ρ * I * S/N – μ * I
𝑑𝑟/𝑑𝑡= μ * I

编写模型内核函数。


> calc4 <- function(rio = 0.3,                 # 接触率,传染率系数
+                   mu=0.1,                    # 治愈率
+                   times = 0:150,             # 病情发展时间
+                   S=1000*1000*10,            # 易感染者1000万人
+                   I=10,                      # 已感染者10人
+                   R=5) {                     # 已移出者5人
+   sir_equations<-function(time,vars,params){
+     with(as.list(c(vars,params)),{
+       N<-S+I+R
+       dS<- -rio*I*S/N
+       dI<-rio*I*S/N - mu*I
+       dR<-mu*I
+       return(list(c(dS,dI,dR)))
+     })
+   }
+   ode(
+     func=sir_equations,
+     y=c(S=S,I=I,R=R),
+     times=times,
+     parms=c(rio=rio,mu=mu)
+   ) %>% as.data.frame()
+ }

当日接触率ρ(rio) = 0.25,日治愈率 μ(mu) =0.1,易感染者S=1000*1000*10,已感染者 I=10,移出者R=100


> sir1<-calc4(rio=0.25,mu=0.1,times=0:150,S=1000*1000*10,I=10,R=100)
> head(sir1)
  time        S        I        R
1    0 10000000 10.00000 100.0000
2    1  9999997 11.61831 101.0789
3    2  9999994 13.49852 102.3324
4    3  9999991 15.68299 103.7887
5    4  9999986 18.22099 105.4808
6    5  9999981 21.16971 107.4466
> sdf4<-melt(sir1,id.vars = c("time"))
> head(sdf4)
  time variable    value
1    0        S 10000000
2    1        S  9999997
3    2        S  9999994
4    3        S  9999991
5    4        S  9999986
6    5        S  9999981
> names(sdf4)<-c("time","type","value")
> draw(sdf4,title="SIR模型",ylab="感染人数",xlab="时间周期")

通过对传播病动力学的基本模型学习,让我了解了传染病的传播基本知识,使我更加坚信坚持“动态清零”的意义是重大的。文本以手写R代码实现模型为主,下一篇文章:用专业工具EpiModel玩转传染病模型

参考文章:

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

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

打赏作者

2022 微软Build After Party:用R语言解读传染病模型

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

关于作者

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

转载请注明出处:
http://blog.fens.me/meeting-ms-build-20220618

前言

近期疫情的反复不断,让北京和上海2个超大城市也受到疫情影响。居家办公,远程开会,已经成为了生活的常态。在这样的背景下,我学习了一些传染病的知识,理解国家坚持动态清零政策的布局,是非常有必要的。

使用R语言结合传染病领域前辈们的专业经验,可以让我们快速上手跨学科的领域,并模拟疫情传播的场景。通过本次分享让大家能感受到疫情传播的可怕,以及我们需要积极的面对,和科学的预防。

目录

  1. 我分享的主题:用R语言解读传染病模型
  2. 会议体验和照片分享

1. 我分享的主题:用R语言解读传染病模型

本次分享的主题,用R语言解读传染病模型。新冠疫情几次变异,极大地影响着我们的正常生活和工作。特别是2022年2月以来的Delta变异株感染,在上海和北京这种人口超大型城市中,有着超强的传染力。政府防疫工作的强力介入,隔离和居家已经是常态了,有新闻指出Delta变异株感染1人可传9人。

在流行病学领域,有几种不同传染病的传播模型,可以模拟病毒的传播过程。本次分享将使用R语言,来给大家演示病毒传播的过程。了解了病毒传播的逻辑,能让我们更加坚定战胜病毒的决心。本次分享的PPT和代码,我上传到了github:https://github.com/bsspirit/infect

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

  1. 传染病模型原理:自由增长模型、SI模型、SIS模型、SIR模型
  2. 用R语言手动实现
  3. 基于EpiModel包的自动化实现
  4. 如何获取新冠数据
  5. 北京的数据带入模型预测

在传染病领域,有4种最基本的传染病模型,分别是自由增长模型、SI模型、SIS模型、SIR模型,这4个模型,分别涉及到现实从得病到治愈再到得病等的病人的状态,通过状态转移人数在计算传播效率。这4个模型,都是可以通过微分方程进行求解的,所以我们可以手动撸代码来计算。具体使用可参考:用R语言解读传染病模型

当然在R语言中,传染病领域专家也提供了,专门的工具包来帮助我们解决传染病的计算和模型的问题,这就是EpiModel包。 EpiModel,提供了用于模拟和分析传染病动力学数学模型的工具,支持的流行病模型类包括确定性隔间模型、随机个体接触模型和随机网络模型。疾病类型包括有和没有人口统计的 SI、SIR 和 SIS 流行病,具有可用于扩展的实用程序,以构建和模拟任意复杂性的流行病模型。 网络模型类基于在 R 的 Statnet 软件套件中实现的时间指数随机图模型 (ERGM) 的统计框架。具体使用可参考:专业工具EpiModel解读传染病模型

2. 会议体验和照片分享

本次活动是微软直通车的活动,主要由微软MVP给大家进行一些技术分享,我们不讲虚的都是干活,边讲PPT,边撸代码。

本次会议报名页: https://mp.weixin.qq.com/s/AXbQO718ZwfYhJ–6bGKyA

2.1 会议主题

MVP嘉宾代表团:由 3位MVP组成,郝冠军,卿毅,张丹。

张丹,用R语言进行量化文本分析,像结构化数据一样来管理文本PPT下载

新冠疫情几次变异,极大地影响着我们的正常生活和工作。特别是2022年2月以来的 Delta 变异株感染,在上海和北京这种人口超大型城市中,有着超强的传染力。政府防疫工作的强力介入,隔离和居家已经是常态了,有新闻指出 Delta 变异株感染1人可传9人。

在流行病学领域,有几种不同传染病的传播模型,可以模拟病毒的传播过程。本次分享将使用R语言,来给大家演示病毒传播的过程。了解了病毒传播的逻辑,能让我们更加坚定战胜病毒的决心。

卿毅,Dynamics 365 与 Power Platform 的集成

低代码/零代码是现在企业数字化转型非常重要的一环。在 Microsoft Build 2022 上,微软发布了 Microsoft 智能数据平台,开发者可以通过低代码连接 Microsoft 智能数据平台的上的数据场景 ,低代码/零代码时代不是取代开发人员 ,而是让通过低代码的方式和 Dynamics 365 业务系统一起创建针对混合办公的协作应用。

郝冠军,.NET 与 Visual Studio 的最新更新

来自产品经理和一线开发人员的一手信息。长年奋战在一线的开发者。目前关注于前端和微服务领域。Angular 技术热爱者,《ASP.NET 本质论》作者。

2.2 相关照片

霸姐,微软MVP大中华区项目负责人。

康康,微软MVP项目助理。

感觉线上的分享还是没有线下分享体验好,似乎少一些沟通的氛围。最后,整个分享结束,感谢组织者刘力科,感谢霸姐支持也都辛苦啦。

转载请注明出处:
http://blog.fens.me/meeting-ms-build-20220618

打赏作者

2021 微软Ignite Post Watching Part:用R语言进行量化文本分析

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

关于作者

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

转载请注明出处:
http://blog.fens.me/meeting-ms-ignite-20211214

前言

疫情反复不断,线下开会越来越难,本次以上线进行分享,参与微软技术直通车。希望疫情早日过去,让大家恢复以往的能见面的交流。

本次的主题是用R语言进行量化文本分析,R语言在被大家广泛的使用,全球社区也在不断的壮大,每天都有好的包被创作出来。本文分析已经是老生长谈的话题,让文本分析过程如结构化数据分析一样的还是比较创新的体验。

目录

  1. 我分享的主题:用R语言进行量化文本分析
  2. 会议体验和照片分享

1. 我分享的主题:用R语言进行量化文本分析

本次分享的主题,是用R语言进行文本分析。在互联网的今天,我们每天都会生产和消费大量的文本信息,如报告、文档、新闻、聊天、图书、小说、语音转化的文字等。海量的文本信息,不仅提供扩宽的研究对象和研究领域,也为商业使用带来了巨大的机会。

R在做数据分析有自己独特的优势,方便、简单、高效。我从2个方面来介绍文本分析,先是jiebaR包进行中文文本分词,然后是quanteda包的进行量化文本分析。本次分享的PPT和代码,我上传到了github:https://github.com/bsspirit/quanteda

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

  1. jiebaR,中文文本分词
  2. quanteda,量化文本分析

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

quanteda,以一种新的方式用结构化数据的方式来管理文本。提出以语料库的形式管理文本,语料库被定义为文本的集合,其中包括特定每个文本的文档级变量,和整个集合的元数据。用户可以轻松地按单词、段落、句子甚至用户提供的分隔符分割文本和标签,按文档级变量将它们分组为更大的文档,形成基于逻辑条件的变量组合。具体使用可参考:http://blog.fens.me/r-word-quanteda/

2. 会议体验和照片分享

本次活动是微软直通车的活动,主要由微软MVP给大家进行一些技术分享,我们不讲虚的都是干活,边讲PPT,边撸代码。

本次活动的官方报名页面:https://www.huodongxing.com/event/2626606217622

2.1 会议主题

MVP嘉宾代表团:由 3位MVP组成,张丹,谢佳标,郝冠军。

张丹,用R语言进行量化文本分析,像结构化数据一样来管理文本PPT下载
在互联网的今天,我们每天都会生产和消费大量的文本信息,如报告、文档、新闻、聊天、图书、小说、语音转化的文字等。海量的文本信息,不仅提供扩宽的研究对象和研究领域,也为商业使用带来了巨大的机会。

量化文本分析(Quantitative Analysis of Textual Data),一种新的方式,用结构化数据的方式来管理文本。quanteda包,提出以语料库的形式管理文本,语料库被定义为文本的集合,其中包括特定每个文本的文档级变量,和整个集合的元数据。用户可以轻松地按单词、段落、句子甚至用户提供的分隔符分割文本和标签,按文档级变量将它们分组为更大的文档,形成基于逻辑条件的变量组合。

谢佳标,《Keras深度学习:入门、实战及进阶》PPT下载

Keras是一个对小白用户非常友好而简单的深度学习框架,它是TensorFlow高级集成API,其特点是能够快速实现模型的搭建,是高效地进行科学研究的关键。本主题将介绍如何进行图像及文本数据预处理,并介绍深度学习常用的DNN、CNN、RNN、GAN等模型原理及Keras案例实现。

郝冠军,在 .NET 6 中应用 OpenTelemetryPPT下载

可观察性是微服务化应用的几个核心特性,OpenTelemetry 延续了 OpenTracing 和 OpenCensus 的发展,成为 CNCF 的针对可观察性的新标准。
该分享将首先介绍 OpenTelemetry 的核心概念,微软作为 OpenTelemetry 的核心成员,.NET 平台对于 OpenTelemetry 提供了优异的支持,这里将基于 .NET 6 介绍如何应用 OpenTelemetry 到 .NET 项目中。

2.2 相关照片

张丹, R语言实践者,北京青萌数海科技有限公司CTO,微软MVP。

10年以上互联网应用架构经验,在R、Java、NodeJS、大数据、数据挖掘等方面有深厚的积累。精通量化投资交易策略,熟悉中国金融二级市场、交易规则和投研体系。 熟悉数据学科方法论,在外汇、海关、区块链等领域均有落地的应用。著有《R的极客理想:量化投资篇》、《R的极客理想:工具篇》、《R的极客理想:高级开发篇》,英文版图书被CRC出版集团引进,在美国发行。个人博客:http://fens.me 。

谢佳标,数据挖掘专家,资深AI技术专家和数据挖掘专家,拥有超过14年的技术研发和管理经验,在数据挖掘和人工智能领域有非常丰富的积累。连续6年(2017-2022)被微软评为最具价值专家(MVP),中国现场统计研究会大数据统计分会首届理事。

郝冠军,郝冠军 10 年微软最有价值专家,多年耕耘在开发前沿,《ASP.NET 本质论》作者,《精通 ASP.NET Core MVC》译者。

感觉线上的分享还是没有线下分享体验好,似乎少一些沟通的氛围。最后,整个分享结束,感谢组织者刘力科,感谢霸姐支持也都辛苦啦。

转载请注明出处:
http://blog.fens.me/meeting-ms-ignite-20211214

打赏作者