美国确诊破十万：这篇文章10天前使用R语言准确预测
本文为作者钱新涵使用R语言进行的新冠病毒预测，原文时间是3月18日，采用的数据集是3月16日（中国数据3月17），文中预测是10天后，恰好是今天。所以我对文章做了翻译，英文原文见（https://towardsdatascience.com/visualize-the-pandemic-with-r-covid-19-c3443de3b4e4）。
以下为译文（3月18发布）：
遵循CDC的建议，我们有两件事情要做：阅读有关COVID-19的新闻，以及对确诊数量增加而感到不知所措。在过去的几周里，情况有多糟？我的手机不断震动着来自世界各地的新闻：首先是从我的家乡中国，然后是亚洲其他地区，欧洲，再到美国。纽约有超过一半的城市已经停止通勤。我们已经用拳头相碰代替握手，现在则更进一步：社交隔离。五天前的3月12日，世界卫生组织宣布COVID-19大流行。迄今为止，已确认来自全球142个国家的160,000多人患有这种疾病。
我觉得有责任使用R创建一个追踪器以解决这一大流行。
我正在使用“ nCOV2019 ”，这是由南方医科大学的Yu Guangchuang Yu博士开发的R软件包。该软件包使我们可以访问所有国家/地区的病例的最新数据和历史数据，在地图上绘制数据并创建各种图形。如果您是像我这样有抱负的数据科学家，请随时安装软件包并按照以下步骤创建视觉效果：
部署程序包
探索数据
创建视觉效果
折线图（按国家/地区）
gif的全球COVID-19增长
影响分析
在接下来的10天里，我们将有多少宗案件？
冠状病毒会影响好莱坞吗？
我们还在吃饭吗？（带有OpenTable数据）
接下来本文就按照上面的步骤进行逐步介绍。
部署程序包
提取此程序包中嵌入的数据的基本功能是：
get_nCov2019()查询在线最新信息load_nCov2019()获取历史数据summary并[访问数据plot在地图上显示数据
安装和部署软件包
remotes::install_github(“GuangchuangYu/nCov2019”)
require(nCov2019)
require(dplyr)
数据概览
x <- get_nCov2019()
y <- load_nCov2019()> xChina (total confirmed cases): 81134last update: 2020–03–17 21:19:04> ynCov2019 historical datalast update: 2020–03–16
保持最新非常重要。简单地打印x和y将刷新数据。
> x['global',]
name confirm suspect dead deadRate showRate heal healRate showHeal1 China 81134 128 3231 3.98 FALSE 68800 84.80 TRUE2 Italy 27980 0 2158 7.71 FALSE 2749 9.82 FALSE3 Iran 16169 0 988 6.11 FALSE 5389 33.33 FALSE4 Spain 11178 0 491 4.39 FALSE 571 5.11 FALSE5 South Korea 8320 0 83 1 FALSE 1401 16.84 FALSE6 Germany 7272 0 17 0.23 FALSE 135 1.86 FALSE7 France 6650 0 148 2.23 FALSE 28 0.42 FALSE8 United States 4687 0 93 1.98 FALSE 74 1.58 FALSE9 Switzerland 2269 0 19 0.84 FALSE 4 0.18 FALSE10 United Kingdom 1950 0 56 2.87 FALSE 52 2.67 FALSE
如何创建数据概览？x[‘global’,]返回最新的全局数据，并根据确认的病例数自动排序。
探索数据
首先，让我们探讨当前数据的整体结构。DataExplorer是一个R软件包，可以快速构建可视化文件。
#explore package
library(DataExplorer)
plot_str(x)
该get_nCov2019()功能获得的数据包括3个列表和5个数据帧。这些是中国乃至世界范围内确诊患者，死亡和康复病例的最新数据。
plot_str(y)
为了研究趋势，我专注于历史数据。通过功能获得的历史数据load_nCov2019()是一个包含3个数据帧的列表。第一个名为“数据”的数据框是城市一级的中国历史数据，包括确诊病例，死亡，康复和嫌疑犯的数量。第二个名为“省”的是省一级的汇总数据。第三个“全球”数据框包括世界各地国家的确诊病例，死亡人数和康复情况。这些历史数据涵盖了从2月15日到最新更新的时间范围。
> summary(x['global',])
name confirm suspect dead deadRate
Length:131 Min. : 1.0 Min. : 0.0000 Min. : 0.00 Length:131
Class :character 1st Qu.: 4.5 1st Qu.: 0.0000 1st Qu.: 0.00 Class :character
Mode :character Median : 30.0 Median : 0.0000 Median : 0.00 Mode :character
Mean : 1223.3 Mean : 0.8626 Mean : 45.44
3rd Qu.: 137.5 3rd Qu.: 0.0000 3rd Qu.: 1.00
Max. :81062.0 Max. :113.0000 Max. :3204.00
showRate heal healRate showHeal
Length:131 Min. : 0.0 Length:131 Length:131
Class :character 1st Qu.: 0.0 Class :character Class :character
Mode :character Median : 0.0 Mode :character Mode :character
Mean : 582.3
3rd Qu.: 3.5
Max. :67023.0
我使用 summary()函数来获取数据的统计概览。在全球范围内，确诊病例的中位数为30，而平均值为1223.3。原因是，离群的感染病例对平均确诊病例有很大影响。
创建视觉效果
可视化是探索性数据分析中的另一种关键方法。借助视觉效果，我们可以轻松地在国家和国际层面观察爆发的趋势。
折线图（按国家/地区）
上面已经提到了如何使用x ['global'，]来获得病例数最多的前十个国家。我们还可以绘制折线图，以查看每个国家/地区的案件随时间增长的情况。
#obtain top 10 country
d <- y[‘global’] #extract global data
d <- d[d$country != ‘China’,] #exclude China
n <- d %>% filter(time == time(y)) %>%
top_n(10, cum_confirm) %>%
arrange(desc(cum_confirm))#plot top 10
require(ggplot2)
require(ggrepel)
ggplot(filter(d, country %in% n$country, d$time > ‘2020–02–15’),
aes(time, cum_confirm, color=country)) +
geom_line() +
geom_text_repel(aes(label=country),
function(d) d[d$time == time(y),]) +
theme_minimal(base_size=14) +
theme(legend.position = “none”)
该图显示了中国境外确诊病例最多的前十个国家。意大利和伊朗是感染最严重的国家，并呈指数增长。同时，韩国通过有效的遏制战略，拉平了曲线并放慢了增长速度。其他一些欧洲国家和美国也看到了成千上万的新病例。
全球COVID-19增长
全球的整体情况如何？我们可以使用plot（）函数轻松绘制全局COVID-19确认的数字地图.
x <-get_nCov2019（）
x
plot（x）＃绘制全局图
那随着时间的变化呢？如果我们可以将更改显示为gif，那将是一个好主意。我绘制了从02-12到03-15的每一天的全局地图，并使用R包magick创建了一个gif。
#visualize global growth over time
library(magick)y <- load_nCov2019()d <- c(paste0(“2020–02-”, 12:29), paste0(“2020–03–0”, 1:9), paste0(“2020–03–1”, 0:5))
img <- image_graph(1200, 700, res = 96)
out <- lapply(d, function(date){
p <- plot(y, date=date,
label=FALSE, continuous_scale=TRUE)
print(p)
})
dev.off()animation <- image_animate(img, fps = 2)
print(animation)
现在我们可以看到冠状病毒是如何在中国开始传播的，并且仅在一个月内就传播到了世界上大多数国家。
影响分析
在接下来的10天里，我们将有多少感染？
接下来，我们可以利用现有数据预测案件的未来增长。
以美国为例。
通过查看图表，我们可以轻松地发现案例呈指数级增长。在这种情况下，我们可以应用对数线性回归来建模和预测增长。由于以下原因，我只考虑100例以上已确诊病例的时间段：
处于早期测试中
早期病例大多与旅行有关，而不是社区传播
usdata <- d %>%filter(d$country == ‘United States’ & d$cum_confirm>100) %>%select(time,cum_confirm)library(forecast)case <- ts(usdata[,2], start=1,frequency = 1)fit <- tslm(log(case) ~ trend)fc <- forecast(fit, h=10)
我们的对数线性模型非常适合该数据集。调整后的R平方为0.9974，表明我们的模型可以解释99％的方差。
Residuals:Min 1Q Median 3Q Max-0.095126 -0.036574 -0.000754 0.036861 0.069923Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 4.311641 0.032840 131.29 < 2e-16 ***trend 0.288296 0.004462 64.61 1.92e-14 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.05336 on 10 degrees of freedomMultiple R-squared: 0.9976, Adjusted R-squared: 0.9974F-statistic: 4175 on 1 and 10 DF, p-value: 1.92e-14
现在我们可以回答这个问题，在10天之内会有多少个案例？
我们可以forecast()用来进行接下来5天的预测（h时段用于设置预测时段）。
Forecasts:Point Forecast Lo 80 Hi 80 Lo 95 Hi 9515 5382.300 4936.042 5868.904 4683.598 6185.23516 7148.845 6541.200 7812.938 6198.094 8245.43617 9495.195 8666.463 10403.174 8199.463 10995.68618 12611.649 11479.941 13854.922 10843.598 14667.98119 16750.966 15204.003 18455.327 14336.190 19572.48420 22248.863 20132.776 24587.365 18948.614 26123.91221 29551.247 26655.278 32761.849 25038.879 34876.81022 39250.374 35286.012 43660.130 33079.244 46572.76423 52132.888 46705.412 58191.072 43692.648 62203.55324 69243.620 61813.326 77567.075 57700.726 83095.642
预测结果表明，如果确诊病例继续呈指数增长，则3天之内将翻一番，而10天之内将达到近7万例。
这就是为什么采取预防措施如此重要的原因。经过数周的逐渐增长（最初似乎可以遏制），感染人数将突然增加。医疗能力将不堪重负，医护人员将处于危险之中。
至关重要的是，我们所有人都要开始远离社会，并避免在太晚之前去公共场所放慢速度。
译者注：结果虽然有所出入，不过整个分析过程非常值得参考。R语言做的gif图也非常漂亮，感谢作者钱新涵，谢谢读者。
往期精选
文章好看点这里