使用 ncov 包获取并分析疫情数据
昨天的推文里，我介绍了开发一个最简单的 R 包的工作流程，相信不少同学已经对 R 包的开发流程有所了解了，今天我们就用这个 ncov 包获取疫情数据然后分析分析吧！
首先我们从 GitHub 上安装或本地安装这个包：
# 从 GitHub 上，需要先安装 devtools 包
devtools::install_github('czxa/ncov')
# 本地安装
install.packages('~/Desktop/ncov_0.0.0.9000.tar.gz', repos = NULL, type = "source")
加载这个包：
library(ncov)
初始化一个 ncov 对象，初始化的过程中所有的数据都会准备好存储在 df 中：
df <- ncov$new()
# 可以查看 df 所属的类
class(df)
#> [1] "ncov" "R6"可以看到 df 是一个 ncov 类，而这个类是基于 R6 类进行封装的，下面再看下这个对象中存储的数据，首先是一些字符串：
# 当前时间：
df$times
#> [1] "截至2月8日12时08分"# 确诊人数
df$confirm
#> [1] "34598"# 疑似人数
df$suspect
#> [1] "27657"# 治愈人数
df$cure
#> [1] "2052"plot() 函数
我在 ncov 类里封装了一个 plot 函数，这个函数封装了 hchinamap 包的 hchinamap 函数，因此可以直接使用 plot() 作用 ncov 对象进行绘图：
# 确诊人数的省份分布：
plot(df,
itermName = "确诊人数",
title = paste0("新冠肺炎确诊人数的分布：", df$times),
subtitle = "数据来源：新浪新闻",
theme = "sandsignika")
这样绘图过于粗糙，我们下面还是从数据出发，手工一步步分析疫情数据。
确诊人数的省份分布
prov_distribution <- df$prov_distribution
prov_distribution
#> # A tibble: 34 x 7
#> name ename value susNum deathNum cureNum city
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <list>
#> 1 北京 beijing 315 157 2 34 <df[,6] [16 × 6]>
#> 2 湖北 hubei 24953 0 699 1119 <df[,6] [17 × 6]>
#> 3 广东 guangdong 1075 195 1 98 <df[,6] [20 × 6]>
#> 4 浙江 zhejiang 1048 0 0 127 <df[,6] [11 × 6]>
#> 5 河南 henan 981 0 4 99 <df[,6] [18 × 6]>
#> 6 湖南 hunan 803 0 0 121 <df[,6] [14 × 6]>
#> 7 重庆 chongqing 426 0 2 31 <df[,6] [39 × 6]>
#> 8 安徽 anhui 733 0 0 47 <df[,6] [16 × 6]>
#> 9 四川 sichuan 363 0 1 51 <df[,6] [21 × 6]>
#> 10 山东 shandong 407 0 0 38 <df[,6] [15 × 6]>
#> # … with 24 more rows
绘制地图就很简单了：
library(sf)
library(tidyverse)
read_sf('china.json') %>%
left_join(prov_distribution) %>%
select(name, value, geometry) -> provmapdata
provmapdata
#> Simple feature collection with 36 features and 2 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: -1000 ymin: 1462 xmax: 9881 ymax: 9848
#> epsg (SRID): 3415
#> proj4string: +proj=lcc +lat_1=18 +lat_2=24 +lat_0=21 +lon_0=114 +x_0=500000 +y_0=500000 +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=m +no_defs
#> # A tibble: 36 x 3
#> name value geometry
#> <chr> <dbl> <MULTIPOLYGON [m]>
#> 1 北京 315 (((6909 6426, 6890 6423, 6890 6423, 6890 6423, 68…
#> 2 天津 88 (((6909 6426, 6933 6438, 6937 6447, 6937 6448, 69…
#> 3 河北 195 (((6941 6461, 6926 6473, 6926 6473, 6926 6473, 69…
#> 4 山西 104 (((6275 5549, 6272 5526, 6260 5505, 6268 5472, 62…
#> 5 内蒙古 52 (((8041 7259, 8035 7233, 8037 7213, 8023 7169, 80…
#> 6 辽宁 99 (((7268 6728, 7284 6732, 7279 6749, 7292 6754, 72…
#> 7 吉林 69 (((8438 6683, 8416 6697, 8431 6708, 8439 6728, 84…
#> 8 黑龙江 295 (((7932 7962, 7932 7962, 7919 7966, 7908 7959, 78…
#> 9 上海 281 (((7881 4314, 7876 4311, 7868 4319, 7809 4314, 77…
#> 10 江苏 439 (((7872 4501, 7821 4512, 7766 4530, 7743 4552, 77…
#> # … with 26 more rows
如果我们想为每个省份添加一个地图的话，我们可以先计算每个省的质心:
provmapdata <- cbind(provmapdata, st_centroid(provmapdata))
glimpse(provmapdata)
#> Observations: 36
#> Variables: 6
#> $ name <chr> "北京", "天津", "河北", "山西", "内蒙古", "辽宁", "吉林", "…
#> $ value <dbl> 315, 88, 195, 104, 52, 99, 69, 295, 281, 439…
#> $ name.1 <chr> "北京", "天津", "河北", "山西", "内蒙古", "辽宁", "吉林", "…
#> $ value.1 <dbl> 315, 88, 195, 104, 52, 99, 69, 295, 281, 439…
#> $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((6909 6426, ...,…
#> $ geometry.1 <POINT [m]> POINT (6764.389 6447.804), POINT (6938…
由于湖北省的确诊数量比其它省份实在是高太多了，所以我们还是把 value 这个变量先分段然后再进行绘图：
provmapdata %>%
mutate(
value = cut(
value,
breaks = c(0.0, 0.9, 9.9, 99.9, 999.9, 9999.9, 99999.9),
labels = c("无", "1～9人", "10～99人", "100～999人", "1000～10000人", "> 10000人")
),
value = as.character(value),
value = case_when(
is.na(value) ~ "无",
T ~ value
)
) -> provmapdata
# 绘制地图
ggplot(provmapdata) +
geom_sf(aes(geometry = geometry, fill = value),
color = "white", size = 0.1) +
geom_sf_text(aes(geometry = geometry.1,
label = name),
family = cnfont, size = 3,
color = "gray") +
scale_fill_viridis_d(name = "确诊人数") +
theme_modern_rc(base_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
labs(title = paste0("新冠肺炎的确诊人数：", df$times),
subtitle = paste0("全国总确诊：", df$confirm, " 例；",
"疑似：", df$suspect, " 例；\n",
" 治 愈：",
df$cure, " 例；",
"死亡：", df$dead, " 例。"),
caption = "数据来源：新浪新闻 ｜ 绘制：TidyFriday")
确诊人数的市级分布
由于我提供的地图数据里面直辖市是作为整体的，所以我手动把获取到的市级数据里面的直辖市数据加总了一下。注意，这里我是使用城市名称的前两个字匹配的，这样匹配的成功率高一些：
city_distribution <- df$city_distribution
city_distribution
# 导出为 csv 文件再手动整理下：
# city_distribution %>%
# select(name, conNum) %>%
# write_csv('city_distribution.csv')
city_distribution <- read_csv('city_distribution.csv') %>%
mutate(name = str_sub(name, 1, 2))
# 市级地图
citymap <- read_sf('china_city.json') %>%
mutate(name = str_sub(name, 1, 2)) %>%
left_join(city_distribution)
# 国界线
cn_boundary <- read_sf('cn_boundary.json')
ggplot() +
geom_sf(data = citymap, aes(fill = log(conNum)),
color = "white", size = 0.05) +
geom_sf(data = cn_boundary, color = "white", size = 0.05) +
scale_fill_viridis_c(name = "确诊人数\n（对数）") +
theme_modern_rc(base_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
labs(title = paste0("新冠肺炎的确诊人数：", df$times),
subtitle = paste0("全国总确诊：", df$confirm, " 例；",
"疑似：", df$suspect, " 例；\n",
" 治 愈：",
df$cure, " 例；",
"死亡：", df$dead, " 例。"),
caption = "数据来源：新浪新闻 ｜ 绘制：TidyFriday") +
coord_sf(crs = "+proj=laea +lat_0=23 +lon_0=113 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
确诊人数的全球分布
这幅图之前画过：
wmap <- read_sf('world.geo.json')
wmap
# 世界各国 ISO 代码对照表
code <- readxl::read_xls('各国国籍代码.xls', skip = 5) %>%
select(id = 代码, name = 内容)
# 合并 wmap 和 othercountry
wdata <- wmap %>%
left_join(code, by = c("id" = "id")) %>%
select(name = name.y) %>%
left_join(df$othercountry) %>%
mutate(
value = case_when(
name == "中国" ~ df$confirm,
T ~ value
),
value = as.numeric(value),
value = cut(
value,
breaks = c(0.0, 0.99, 9.9, 49.9, 99.9, 999.9, 100000),
labels = c("无", "1～9人", "10～49人", "50～99人", "100～1000人", "> 1000人")
) %>% as.character()
) %>%
mutate(
value = case_when(
is.na(value) ~ "无",
T ~ value
),
value = factor(
value,
levels = c("无", "1～9人", "10～49人", "50～99人", "100～1000人", "> 1000人")
)
)
# 绘制世界地图
ggplot() +
geom_sf(data = wdata,
aes(fill = value),
size = 0.1, color = "white") +
theme_modern_rc(base_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont) +
coord_sf(crs = "+proj=eck4") +
worldtilegrid::theme_enhance_wtg() +
scale_fill_viridis_d(name = "确诊人数",
direction = 1) +
labs(title = "世界各国感染新型冠状病毒肺炎确诊人数分布",
subtitle = df$times,
caption = "数据来源：新浪新闻 | 绘制：TidyFriday")
历史数据趋势
突然发现我忘记把历史数据封装进 ncov 对象里面了！没事，还可以使用 jsondata 进行提取：
library(lubridate)
library(scales)
df$jsondata$data$historylist %>%
as_tibble() %>%
type_convert() %>%
select(
日期 = date,
确诊 = cn_conNum,
疑似 = cn_susNum,
治愈 = cn_cureNum,
死亡 = cn_deathNum
) %>%
gather(确诊, 疑似, 治愈, 死亡,
key = "kind", value = "value") %>%
mutate(日期 = ymd(paste0("2020.", 日期))) %>%
ggplot(aes(日期, value, color = kind)) +
geom_line() +
geom_point() +
scale_color_viridis_d(option = "A", begin = 0.3, end = 1) +
theme_modern_rc(base_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont) +
labs(x = "", y = "", title = "新型冠状病毒肺炎疫情发展情况",
subtitle = paste0(df$times, "，确诊人数：",
df$confirm, "人"),
caption = "数据来源：新浪新闻 | 绘制：TidyFriday",
color = "") +
scale_x_date(breaks = date_breaks('3 days'),
labels = date_format("%m-%d"))
新冠病毒的传染能力：R0
前面已经介绍过 R0 的计算了，这里就不再介绍了。我们来使用最新的数据计算 R0:
library(deSolve)
R0 <- function(days){
historylist <- suppressMessages(
df$jsondata$data$historylist %>%
type_convert()
)
# 感染人数的序列
da <- historylist %>% pull(cn_conNum) %>% rev()
# 选择从第六天开始的：
Infected <- da[6:days]
# 天数
Day <- 1:(length(Infected))
# 中国的总人口
N <- 1400000000
SIR <- function(time, state, parameters) {
par <- as.list(c(state, parameters))
with(par, {
dS <- -beta/N * I * S
dI <- beta/N * I * S - gamma * I
dR <- gamma * I
list(c(dS, dI, dR))
})
}
Day = 1:(length(Infected))
init = c(S = N-Infected[1], I = Infected[1], R = 0)
RSS = function(parameters) {
names(parameters) <- c("beta", "gamma")
out <- ode(y = init, times = Day, func = SIR, parms = parameters)
fit <- out[ , 3]
sum((Infected - fit)^2)
}
Opt = optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1))
Opt_par = setNames(Opt$par, c("beta", "gamma"))
return(Opt_par["beta"] / Opt_par["gamma"])
}
library(purrr)
R0DF <- tibble(days = 7:length(da)) %>%
mutate(
R0 = map_dbl(days, R0),
date = ymd("2020-01-11") + days(days)
)
print(R0DF, n = Inf)
#> # A tibble: 22 x 3
#> days R0 date
#> <int> <dbl> <date>
#> 1 7 1.94 2020-01-18
#> 2 8 2.85 2020-01-19
#> 3 9 2.94 2020-01-20
#> 4 10 2.63 2020-01-21
#> 5 11 2.67 2020-01-22
#> 6 12 2.53 2020-01-23
#> 7 13 2.46 2020-01-24
#> 8 14 2.45 2020-01-25
#> 9 15 2.45 2020-01-26
#> 10 16 2.41 2020-01-27
#> 11 17 2.43 2020-01-28
#> 12 18 2.40 2020-01-29
#> 13 19 2.34 2020-01-30
#> 14 20 2.29 2020-01-31
#> 15 21 2.23 2020-02-01
#> 16 22 2.17 2020-02-02
#> 17 23 2.12 2020-02-03
#> 18 24 2.07 2020-02-04
#> 19 25 2.03 2020-02-05
#> 20 26 1.99 2020-02-06
#> 21 27 1.95 2020-02-07
#> 22 28 1.91 2020-02-08
从数据上就能看到，R0 下降的很快，而且很稳定，这是因为我使用的是增量数据（并没有固定窗口）。再绘图展示：
library(latex2exp)
ggplot(R0DF, aes(x = date, y = R0)) +
geom_point(color = "#18BC9C") +
geom_line(color = "#18BC9C") +
labs(x = "", y = "基本传染数",
title = TeX("新型冠状病毒的基本传染数：R_0"),
subtitle = TeX("R_0：表示平均有多少个健康的人被患者感染"),
caption = "TidyFriday Project") +
scale_x_date(breaks = date_breaks('3 days'),
labels = date_format("%m-%d")) +
theme_modern_rc(base_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont)
预测
一般来说，如果一个传染病的 R0 低于 1，那么这个传染病就会自动消失；如果 R0 为 0 就说明这个病没有传染能力。所以我们关注两个时点，第一，什么时候 R0 会低于 1；第二，什么时候 R0 会达到 0。0 之后的预测就没有意义了。
sweep 包提供了一套时间序列预测的整洁工具（tidy tool），首先我们加载需要的一些包：
library(sweep)
library(forecast)
library(timetk)
library(sweep)
去除 R0DF 的第一列：
R0DF <- R0DF[,2:3]
建模预测：
# 根据 R0DF 生成一个时间序列，start 不重要，等下再处理
R0TS <- tk_ts(R0DF, start = 2020, freq = 1, silent = TRUE)
# 使用 ETS (指数平滑模型) 模型：
fit_ets <- R0TS %>%
ets()
# 预测未来 50 天的：
fcast_ets <- fit_ets %>%
forecast(h = 50)
绘图展示预测结果：
sw_sweep(fcast_ets) %>%
mutate(index = index - 2020) %>%
mutate(date = ymd("2020-01-18") + days(index)) %>%
ggplot(aes(x = date, y = R0)) +
geom_ribbon(aes(ymin = lo.95, ymax = hi.95), fill = "#D5DBFF", color = NA, size = 0) +
geom_ribbon(aes(ymin = lo.80, ymax = hi.80), fill = "#596DD5", color = NA, size = 0, alpha = 0.8) +
geom_line(size = 1) +
geom_hline(yintercept = 0,
color = "#18BC9C", linetype = 2) +
geom_hline(yintercept = 1,
color = "#E31A1C", linetype = 2) +
labs(
title = "新型冠状病毒 R0 的变化和预测",
x = "", y = TeX("R_0"),
subtitle = TeX("R_0：表示平均有多少个健康的人被患者感染"),
caption = "TidyFriday Project") +
scale_color_tq() +
scale_fill_tq() +
scale_x_date(breaks = date_breaks('7 days'),
labels = date_format("%m-%d")) +
theme_modern_rc(base_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont)
❝附件链接：https://t.zsxq.com/iMfEIE6 ❞
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2020-02-09，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录
相关产品与服务
对象存储
对象存储（Cloud Object Storage，COS）是由腾讯云推出的无目录层次结构、无数据格式限制，可容纳海量数据且支持 HTTP/HTTPS 协议访问的分布式存储服务。腾讯云 COS 的存储桶空间无容量上限，无需分区管理，适用于 CDN 数据分发、数据万象处理或大数据计算与分析的数据湖等多种场景。