新型冠状病毒的传染性如何?

新型冠状病毒的传染性如何?

昨天看到一篇博客,关于计算新型冠状病毒的 R0(基本传染数) 的,感觉很有意思,于是就想复现一遍,原文在这里:Epidemiology: How contagious is Novel Coronavirus (2019-nCoV)?

注意:本文中所有的计算结果仅供演示,请勿作为专业的研究结果。

从新浪新闻提供的数据接口爬取数据

新浪新闻提供的接口在这里:https://interface.sina.cn/news/wap/fymap2020_data.d.json 是个 json 格式的数据,我们可以使用 jsonlite 包处理:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(jsonlite)
library(tidyverse)
jsondata <- fromJSON('https://interface.sina.cn/news/wap/fymap2020_data.d.json')
# 当前时间
(times <- jsondata$data$times)
#> "截至2月5日10时30分"

# 确诊数量
(confirm <- jsondata$data$gntotal)
#> [1] "24363"

# 死亡数量
(dead <- jsondata$data$deathtotal)
#> [1] "490"

# 疑似数量
(suspect <- jsondata$data$sustotal)
#> [1] "23260"

# 治愈数量
(cure <- jsondata$data$curetotal)
#> [1] "892"

新型冠状病毒肺炎的省份分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
prov_distribution <- jsondata$data$list %>%
as_tibble()
prov_distribution
#> # A tibble: 34 x 7
#> name ename value susNum deathNum cureNum city
#> <chr> <chr> <chr> <chr> <chr> <chr> <list>
#> 1 北京 beijing 228 0 1 23 <df[,6] [15 × 6]>
#> 2 湖北 hubei 16678 0 479 520 <df[,6] [17 × 6]>
#> 3 广东 guangdong 870 136 0 32 <df[,6] [20 × 6]>
#> 4 浙江 zhejiang 895 0 0 70 <df[,6] [11 × 6]>
#> 5 河南 henan 764 0 2 41 <df[,6] [18 × 6]>
#> 6 湖南 hunan 661 0 0 35 <df[,6] [14 × 6]>
#> 7 重庆 chongqing 366 0 2 14 <df[,6] [37 × 6]>
#> 8 安徽 anhui 530 0 0 20 <df[,6] [16 × 6]>
#> 9 四川 sichuan 301 0 1 23 <df[,6] [21 × 6]>
#> 10 山东 shandong 298 0 0 13 <df[,6] [15 × 6]>
#> # … with 24 more rows

前面的地图画的太多了,这里就不画地图了🥱。

新型冠状病毒肺炎的城市分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
city_distribution <- prov_distribution %>%
select(city) %>%
unnest(city) %>%
type_convert()
#> # A tibble: 404 x 6
#> name conNum susNum cureNum deathNum mapName
#> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 海淀区 42 0 0 0 海淀区
#> 2 怀柔区 3 0 0 0 怀柔区
#> 3 延庆区 1 0 0 0 延庆区
#> 4 丰台区 18 0 3 0 丰台区
#> 5 大兴区 28 0 0 0 大兴区
#> 6 东城区 2 0 0 0 东城区
#> 7 昌平区 13 0 0 0 昌平区
#> 8 西城区 28 0 0 0 西城区
#> 9 朝阳区 26 0 0 0 朝阳区
#> 10 石景山区 6 0 0 0 石景山区
#> # … with 394 more rows

同样我也不再画市级的地图了,大家可以自行尝试下。

新型冠状病毒肺炎在国外的分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
othercountry <- jsondata$data$otherlist %>%
as_tibble()
#> # A tibble: 31 x 5
#> name value susNum deathNum cureNum
#> <chr> <chr> <chr> <chr> <chr>
#> 1 德国 12 3 0 0
#> 2 比利时 1 0 0 0
#> 3 西班牙 1 0 0 0
#> 4 俄罗斯 2 0 0 0
#> 5 柬埔寨 1 0 0 0
#> 6 印度 3 0 0 0
#> 7 阿联酋 5 0 0 0
#> 8 斯里兰卡 1 0 0 0
#> 9 泰国 25 0 0 8
#> 10 韩国 18 0 0 1
#> # … with 21 more rows

我们再绘制一幅世界地图

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
# 世界地图
library(sf)
library(ggplot2)
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(othercountry) %>%
mutate(
value = case_when(
name == "中国" ~ 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 = times,
caption = "数据来源:新浪新闻 | 绘制:TidyFriday")

历史数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
library(lubridate)
historylist <- jsondata$data$historylist %>%
mutate(
date = ymd(paste0("2020.", date))
) %>%
as_tibble() %>%
type_convert()
historylist

#> # A tibble: 25 x 10
#> date cn_conNum cn_deathNum cn_cureNum cn_susNum wjw_susNum
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2020-02-04 24363 490 892 23260 3971
#> 2 2020-02-03 20471 425 632 23214 5072
#> 3 2020-02-02 17238 361 475 21558 5173
#> 4 2020-02-01 14411 304 328 19544 4562
#> 5 2020-01-31 11821 259 243 17988 5019
#> 6 2020-01-30 9720 213 171 15238 4812
#> 7 2020-01-29 7736 170 124 12167 4148
#> 8 2020-01-28 5974 132 103 9239 3248
#> 9 2020-01-27 4515 106 60 6973 2077
#> 10 2020-01-26 2744 80 51 5794 3806
#> # … with 15 more rows, and 4 more variables: wuhan_conNum <dbl>,
#> # wuhan_deathNum <dbl>, wuhan_cureNum <dbl>, wuhan_susNum <lgl>

可以绘制一幅折线图展示这一数据:

绘制代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(scales)
historylist %>%
select(
日期 = date,
确诊 = cn_conNum,
疑似 = cn_susNum,
治愈 = cn_cureNum,
死亡 = cn_deathNum
) %>%
gather(确诊, 疑似, 治愈, 死亡,
key = "kind", value = "value") %>%
ggplot(aes(日期, value, color = kind)) +
geom_line() +
geom_point() +
awtools::a_secondary_color() +
theme_ipsum(base_family = cnfont) +
labs(x = "", y = "", title = "新型冠状病毒肺炎疫情发展情况",
subtitle = paste0(times, ",确诊人数:", confirm, "人"),
caption = "数据来源:新浪新闻 | 绘制:TidyFriday",
color = "") +
scale_x_date(breaks = date_breaks('3 days'),
labels = date_format("%m-%d"))

我们可以把确诊人数单独拿出来观察:

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
historylist %>%
select(
日期 = date,
确诊 = cn_conNum
) %>%
ggplot(aes(日期, 确诊)) +
geom_point(color = "black") +
geom_smooth(method = 'lm') +
theme_ipsum(base_family = cnfont) +
labs(x = "", y = "", title = "新型冠状病毒肺炎疫情发展情况",
subtitle = paste0(times, ",确诊人数:", confirm, "人"),
caption = "数据来源:新浪新闻 | 绘制:TidyFriday") +
scale_y_log10() +
scale_x_date(breaks = date_breaks('3 days'),
labels = date_format("%m-%d"))

SIR 模型

计算最新的 R0

关于 SIR 模型我就不班门弄斧地讲解了,大家有兴趣地可以阅读原文,下面我们根据原文提供的代码计算 R0,为了得到更准确的计算,我们从第六个观测值开始计算(前五个都是 41):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# 感染者人数
da <- historylist %>% pull(cn_conNum) %>% rev()
Infected <- da[6:length(da)]
# 天数
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))
})
}

library(deSolve)
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$message
#> "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

是个收敛的结果,表明最优化成功。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
#> beta gamma
#> 0.6701221 0.3298780

# 时间:1 ~ 70 天
t <- 1:70
fit <- data.frame(ode(y = init, times = t, func = SIR, parms = Opt_par)) %>%
as_tibble()
fit
#> # A tibble: 70 x 4
#> time S I R
#> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1399999959 41 0
#> 2 2 1399999929. 53.8 17.2
#> 3 3 1399999890. 70.5 39.7
#> 4 4 1399999838. 92.5 69.2
#> 5 5 1399999771. 121. 108.
#> 6 6 1399999682. 159. 159.
#> 7 7 1399999566. 209. 225.
#> 8 8 1399999414. 274. 313.
#> 9 9 1399999214. 359. 427.
#> 10 10 1399998952. 471. 577.
#> # … with 60 more rows

这个图的绘制代码为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
fit %>% 
gather(S, I, R, key = "kind", value = "value") %>%
ggplot(aes(time, value, color = kind)) +
geom_point() +
geom_line() +
awtools::a_flat_color(
breaks = c("S", "I", "R"),
labels = c("易感", "感染", "康复")
) +
theme_ipsum(base_family = cnfont) +
scale_y_continuous(
breaks = c(0, 5e+08, 1e+9, 1.5e+9),
labels = c("0", "5 亿", "10 亿", "15 亿"),
expand = c(0, 0.2e+9)
) +
labs(title = "新型冠状病毒肺炎疫情的 SIR 模型",
x = "天数", y = "人数", color = "")

根据模型拟合结果就可以计算 R0 了:

1
2
3
4
R0 <- setNames(Opt_par["beta"] / Opt_par["gamma"], "R0")
R0
#> R0
#> 2.031424

可以看到根据这个模型,使用最新数据计算的 R0 是 2.03。

如果接下来任由病毒的传播,那么最大感染人数:

1
2
3
4
5
fit[fit$I == max(fit$I), c("time", "I"), drop = FALSE]
#> # A tibble: 1 x 2
#> time I
#> <dbl> <dbl>
#> 1 51 221643584.

最大感染人数为 2.2 亿(将在 51 天的时候达到,从 1 月 16 号开始,也就是 3 月 7 号)!如果按照 2% 的死亡率,将会造成 443 万人的丧生!

1
2
max(fit$I) * 0.02
[1] 4432872

观察 R0 的变化趋势

当然这里计算的一切都不可能发生,从数据来看最近的 R0 正在快速的下降,我们可以看一下前几天的数据计算的结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
# 首先写一个计算 R0 的函数:
R0 <- function(days){
Infected = da %>%
.[6:days]
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)
RODF <- tibble(days = 7:length(da)) %>%
mutate(
R0 = map_dbl(days, R0),
date = ymd("2020-01-11") + days(days)
)
RODF

#> # A tibble: 19 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

绘图展示:

绘图代码:

可以使用 latex2exp 包方便地在绘图中使用 LaTeX 公式:

1
2
3
4
5
6
7
8
9
10
library(latex2exp)
ggplot(RODF, aes(x = date, y = R0)) +
geom_point(color = "black") +
geom_line(color = "black") +
labs(x = "", y = "基本传染数",
title = TeX("新型冠状病毒的基本传染数:R_0"),
subtitle = TeX("R_0:表示平均有多少个健康的人被患者感染")) +
theme_ipsum(base_family = cnfont) +
scale_x_date(breaks = date_breaks('3 days'),
labels = date_format("%m-%d"))

本文所需的材料可以前往知识星球下载:https://t.zsxq.com/MVZn2BY

# R

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×