原创：MarTechApe

Contributors | Vivian老师，大妍，小虎 — Contributors | Vivian老师，大妍，小虎

一场发生于12月末的新型冠状病毒肺炎以武汉为中心，在几天的时间内迅速扩散至全国，甚至海外部分地区。如今，距离疫情最初的爆发已经过去了一个多月的时间，在医护人员的努力下，战疫已经有了初步的成果，而疫情相关的新闻报道依旧是大家关注的重点。

这场惊动全球的疫情一个多月来究竟演变如何？之前的专家预测与现在的实况相比是否准确？疫情拐点何时到来？借助公开数据，我们做了一些分析，看完本文，相信你也可以自己动手试试分析一下。

接下来，小编就带你用数据来分析和预测一下这场战疫的情况。

1.数据准备和初印象

一周前，余光创博士在Github上发布了名为nCoV2019的R语言包，可以快速获取本次疫情相关统计数据。这里，我们主要使用R作为分析工具，（R是属于GNU系统的一个自由、免费、源代码开放的软件，主要用于统计分析、绘图的语言和操作环境。所有R的函数和数据集是保存在程序包里面的。只有当一个包被载入时,它的内容才可以被访问。）所以，当我们在使用这个nCoV2019的语言包之前，我们需要在R里面安装remotes package并调用这个包。remotes这个package的作用就是帮大家下载并安装Github上的R package。

在R中，安装package只需要进行一次：

install.packages("remotes")

安装成功后，每次做分析前，需要使用library function来调用package：

library(remotes)

然后在R里面运行以下语句来下载余博士创建的疫情相关的数据

remotes::install_github("GuangchuangYu/nCov2019")

数据下载完成后，我们需要来调用这个包，获取数据进行分析。

library(nCov2019)

x <- get_nCov2019()

我们将这个package中所有和疫情相关的数据信息存放在x这个变量中，如果我们想要快速地查看x中包含了哪些元素，我们可以使用$来调用其中的元素并查看，举个例子：

> x$lastUpdateTime

[1] "2020-02-17 09:18:14"

通过这个语句，我们可以查看到这个R package中最近一次更新数据是在中国时间的2020年的2月17日早上9：18分。

如果查看chinaTotal这个数据列表，我们会发现，其中包括了累计确诊、现存疑似、死亡和治愈人数。同样的道理，可以通过x$chinaAdd来查看当日的新增数据。

截至2020-02-17，上午9时18分，国内累计确诊人数为70635例，现存疑似为7264例，死亡为1772例，治愈为10853例。而2月16日官方公布数据为，累计确诊病例68592例, 累计治愈人数为9671，现存疑似人数为8228。将两天的数据相对比，我们能看出确诊人数增加（2043）的同时，治愈病例也在稳步增加（1182），并且累计疑似病例也出现了下降。这也验证了专家所预测的重要拐点：通过封城等阻断传染源减少了病毒的进一步传播，并且火神山和雷神山的建设成功大大加快了病例确诊的速度，使大批疑似病例得到确诊，开始治疗。这也使得我们所看到的确诊病例数据出现了大幅度提升。（source：丁香园-丁香医生）

以上我们所使用的的方法是一种非常快速地了解数据结构的方法，但是也许你会觉得这并不直观，也无法全面地了解其中元素之间的关联，这里给大家安利R里面一个叫做DataExplorer的package来一目了然地查看package中包含了哪些信息

library(DataExplorer)

plot_str(x）

从以上图表中，我们可以看到这个Package中一共包括2个list（列表）和7个dataframe（数据框）。在R中，List是最宽泛的一种数据结构，它包括的数据类型可以不一致，长度也可以不一致。可以由向量，矩阵，数组，数据框，函数，甚至是列表组成。而Dadaframe数据框里面列与列之间的数据类型可以不一致，但同一列中的数据类型必须一致，且每一列长度必须相同（有相同的行数）。

其中主要的Dataframe是 chinaDayList，我们可以通过names()这个语句看到，这个表中包含了每日确诊数、疑似数、死亡数、治愈数、死亡率与治愈率。

> names(x$chinaDayList)

[1] "confirm" "suspect" "dead" "heal" "deadRate" "healRate" "date" "confirm_cuts" "suspect_cuts" "dead_cuts" "heal_cuts"

我们还可以通过命令x[ ]，x[‘global’, ]调取最新全国甚至是全世界的数据，将数据精确到各国家和各省份。

通过数据列表，我们可以看到中国和其周边国家病例数量最为显著，中国国内以湖北、广东、河南、浙江、湖南等省病情最为严重。

以上是对整个package的一个第一印象，当我们对package中的数据结构和元素有了一个初步了解以后，我们就可以对数据进行一些初步的分析了。

2. 探索性数据分析

知乎中有一篇关与探索性数据分析的文章写道探索性数据分析（Exploratory Data Analysis, 简称EDA）像侦探工作，你不知道自己会找到什么，所以你也不会做过多的假设，通过可视化工具对数据进行一次又一次的检视来找到线索，并对结果保持开放的心态。（source：https://zhuanlan.zhihu.com/p/28470886）

在之前的文章中我们也介绍过， R是由统计学家开发，可视化能力很强的一种分析工具。这非常有利于我们进行一些探索性的数据分析。探索性数据分析最重要的一个目的就是培养你对数据的直觉。既然如此，不妨从一些描述统计开始。我们可以通过summary()来实现。

从这个总结中，我们能在这个函数获取描述性统计量，比如最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计。举个例子，从时间（date）和累计疑似病例（suspect）来看，我们能很容易地理解：从1月13日至2月17日的35天的时间（Length：35），最少的累计疑似病例数为0例，最多的累计疑似病例为28942例，其中平均值约为为10805例，中间值为8969例。这些数字也体现出了全国病情严重程度。描述统计可以帮我们快速的找到数据中的几个关键节点，但是单纯的依赖描述统计可以得到的有效信息有限。

除了描述统计以外，可视化是EDA的精髓，在R中，最常用的进行数据可视化的package是ggplot2, 相信使用过R的小伙伴对这个package一定不陌生。接下来，就让我们用ggplot来画一下每日累计确诊人数的走势。

ggplot(d$chinaDayList,

aes(as.Date(date, "%M.%d"),

as.numeric(confirm))) +

geom_col(fill='steelblue')

这里横轴x为日期，纵轴y为确诊病例数量。以上的柱状图更能清晰且直观地为我们展现出了累计确诊病例的直线式上升趋势。此外，在2月12日我们可以看到一个明显的激增，随后增幅和之前几天的增幅基本一致。当看到一个如此明显的上涨后，我们带着疑问，通过一些调查研究不难得知，激增是由于卫健委将临床诊断病例数纳入确诊病例数导致的。

上文提到，我们可以在这个package中获取最新全国各个省份以及世界其他地区的疫情数据，那么疫情按照地域的分布图我们可以通过以下指令来实现。

首先，我们来观察一下世界地图，这里我们使用的画图工具为“plot”。在R里，plot()函数是一种常用的绘图函数，用其可以绘制散点图、曲线图等。

> require(nCov2019)

> x = get_nCov2019()

> plot(x)

我们再从世界缩小到全国地图，进一步来观察一下全国的病情分布情况。在绘制中国地区疫情分布图时，我们需要借助mapdata这个package来划分省界。

install.packages("mapdata")

library(mapdata)

remotes::install_github("GuangchuangYu/chinamap") #下载地图

下载好中国地图包之后，我们便可以通过将region设为china，得到全国的病情分布情况。

plot(x, region='china', chinamap=cn)

3. 预测分析

疫情的发展牵动着我们每一个人的心，相信大家和小编一样，一直在实时关注着疫情的发展状况，也很想能够拥有一个水晶球，知道疫情在未来将会如何发展？预知这场灾祸什么时候可以结束。下面让我们试着利用现有的数据对疫情进行一些预测。通过查阅资料，发现来自中南大学和其他高校的学者于2月1日发布了一篇学术研究，研究通过1月22日至2月1日的确诊数据建立了数学模型，并对疫情走向进行了预测分析，论文预测确诊数高峰将在2月2日至2月5日出现。

距离论文的发布已经过去了一段时间，现实的数据是否验证了这个结论呢？

我们先用ggplot来绘制一下新增确诊人数的折线图，发现如专家所预测在2月4日出现了确诊数的小高峰，在2月4日后确诊数开始波动下降。但是由于2月12日卫健委将临床诊断病例数纳入确诊病例数导致了确诊人数的又一个高峰。

ggplot(x$dailyNewAddHistory,

aes(as.Date(date, "%m.%d"),

as.numeric(country))) +

geom_line()

由于2月12日之前与之后，确诊病例的标准不同，在进行预测前，我们需要对数据进行一些预处理，使整个时间序列上的数据点遵循相同的确诊标准，具有可比性。根据丁香医生的疫情日报显示，如果确诊标准不变，2月12-13日新增确诊人数分别为1508和1998人。这也对应了根据确诊的标准不同，新增数量出现了剧增的情况。

在遵循相同的确诊标准对数据进行处理后，我们对新增确诊人数进行了重新绘制。图中红色折线代表不含临床诊断的确诊数，灰色折线代表更新标准后的确诊数。我们可以发现，2月4日以后每日增加的确诊数在稳步下降。

ConfirmOld <-cbind(x$chinaDayAddList["date"],x$chinaDayAddList["confirm"])

ConfirmOld$new <- ConfirmOld$confirm

ConfirmOld$new[24:25] <- c(1508,1998)

ggplot(ConfirmOld,

aes(as.Date(date, "%m.%d"))) +

geom_line(aes(y=ConfirmOld$confirm), color = "grey") +

geom_line(aes(y=ConfirmOld$new), color = "red") +

geom_point(aes(y=ConfirmOld$confirm),color = "grey" ) +

geom_point(aes(y=ConfirmOld$new), color = "red" )

对数据进行一个调整之后，我们可以使用R的forecast package对全国每日的新增人数进行指数平滑的预测。当我们的数据是一个没有季节性变动，恒定增加的时间序列时，我们可以使用指数平滑的方法进行短期预测。对于简单的指数平滑，我们可以使用forecast包中的HoltWinters()函数进行预测。

case <- ts(adjusted)

plot.ts(case)

forecast <- HoltWinters(case, beta=FALSE, gamma=FALSE)

结果如下所示，在简单的指数平滑预测中，平滑度由参数alpha控制（alpha取值0-1）。alpha越大，说明预测值越依赖于最新的数据点。在我们的例子中，alpha值高达0.94，说明预测高度依赖最新一天（2/16）的新增确诊人数。