作者：Anderson

翻译：郑金鑫  中国疾病预防控制中心，流行病与卫生统计专业

目前许多国家都出现相应的疫苗问题，丹麦爆发了5例麻疹后，疫苗问题再次成为热点。与日本爆发的麻疹疫情相比，丹麦的疫情不算严重。但是当地卫生部仍然警惕，可能会发生潜在的疾病流行。在这里，我们将使用Shiny创建一个应用程序，描述疫苗接种对传染病的影响过程。在Shiny应用程序中，我们增加了参数选项，用户可自己调试不同参数，以观察疫苗接种对传染病暴发的影响。希望该程序能够加深人们对疫苗接种的认识。

SIR疾病传播类型

susceptible-infectious-recovered （SIR）疾病传播模型是最简单，最经典的一种刻画疾病传播的模型，它可成功用于描述麻疹，流感，水痘，腮腺炎和风疹的疾病爆发。 SIR模型可通过改变参数进行扩展，以适应早期感染或疫苗接种引起的疾病传播。已有大量文献记录了可在R中实现SIR模型，因此在本文我们将在前期工作基础上，基于SIR模型添加一个额外的参数以描述接种疫苗的个体在传播中的作用，并在Shiny应用程序中实现所有内容。

我们将使用4个参数的SIR模型，假设某一群体中的人群可分为4类：易感人群（S），已感染人群（I），恢复人群（R）和免疫接种人群（V）个人。已经感染的人群在开始阶段非常少，因为疾病的传播是由具有传染性的病人传播至他人导致感染人群增加，R中的初始人数为0，因为开始疾病传播开始阶段，并没有发展到治愈阶段。剩下的人群都属于S（之前没有患过疾病）或V（接种疫苗/具有免疫力）。设置如下图所示。将V添加到SIR模型中可以减少疾病的传播，因为如果疾病已经免疫或接种疫苗，它就不会感染与之接触的个体。

SIR模型的这种扩展可用于模拟与解释传染性疾病的传播，在模型中，我们将忽视死于其他原因的疾病，新疫苗接种人数及该人群新增人口也一并忽略。该模型的最重要的前提条件是，总人口规模在一段时间内是恒定的，而且传播率β和γ在这段时间内是恒定的，并且该人群中的任何人都可以与其他任何人接触。要计算爆发的最后效果，我们需要为模型设置一些初始参数。直接影响模型的参数有：

• β - 从S到I的转变率。该比率定义为基本生殖数R0除以感染期（即，每个病人每天有效接触的平均人数）

• γ-从I到R的转变率。该值是疾病期的倒数，因为一旦疾病期结束，该个体将自动转移到R组。

因此，我们在shiny应用程序中，设置了以下参数。

• 基本生殖数R0–每个感染者可能感染的平均个体数，即疾病的传染性

• 感染期

• 人口数量

• 最初感染的人数

• 具有免疫个体的比例V. 如果不是100％，则该百分比应乘以疫苗有效性。

• 要考虑的时间范围

SIR 公式

下面的代码中显示了这些初始参数的示例

 1# Set parameters 2timeframe <- 200    # Look at development over 200 days 3pinf <- 5           # Initial number of infected persons 4popsize <- 5700000  # Population size (5.7 mio in Denmark) 5pvac <- .90         # Proportion vaccinated 6vaceff <- .95       # Effectiveness of vaccine 7connum <- 15        # R0, the reproductive number. 15 is roughly measles 8infper <- 14        # Infection period, 14 days 910# Compute the proportions in each of the compartments at11# the initial outbreak12init <- c(S = 1 - pinf / popsize - pvac * vaceff,13          I = pinf / popsize,14          R = 0,15          V = pvac * vaceff 16         )1718# First set the parameters for beta and gamma19parameters <- c(beta = connum / infper,20                gamma = 1 / infper)21## Time frame22times      <- seq(0, timeframe, by = .2)

基于上面概述的模型，我们现在可以设置以下一组耦合微分方程

如果βS（t）-γ> 0那么受感染的个体数量将增加且疾病将暴发。如果βS（t）-γ<0，则感染的数量将减少，且疾病停止流行。请注意，由于S（t）随着时间的推移而减少，暴发最终将自行消失，因为当前大多数人口集中在恢复组中，从易感人群组进入回复组，不在感染疾病。现在我们准备定义描述传染病传播过程的微分方程。将传染的过程公式化，当前的时间点，time，当前状态 state（即，4个人群）和参数parameters。

 1## Create a SIR function with an extra V component 2sirv <- function(time, state, parameters) { 3  with(as.list(c(state, parameters)), { 4    dS <- -beta * S * I 5    dI <-  beta * S * I - gamma * I 6    dR <-                 gamma * I 7    dV <- 0 8    return(list(c(dS, dI, dR, dV))) 9  })10}

一旦我们有了sirv（）函数，我们就可以使用deSolve包中的ode（）函数来求解耦合微分方程。

 1library("deSolve") 2 3## Solve the set of coupled differential equations 4out <- ode(y = init, 5           times = times, 6           func = sirv, 7           parms = parameters 8          ) 9head(as.data.frame(out))10##   time         S            I            R     V11## 1  0.0 0.1449991 8.771930e-07 0.000000e+00 0.85512## 2  0.2 0.1449991 8.920417e-07 1.263739e-08 0.85513## 3  0.4 0.1449991 9.071419e-07 2.548870e-08 0.85514## 4  0.6 0.1449990 9.224976e-07 3.855756e-08 0.85515## 5  0.8 0.1449990 9.381132e-07 5.184763e-08 0.85516## 6  1.0 0.1449990 9.539932e-07 6.536268e-08 0.855

基本生殖数R0的值取决于传染病类型的种类，下表中看到.R0的大值意味着该疾病具有高度传染性。

常见传染病的基本生殖基数R0

现在具备一切条件后，我们只需要将它整合到Shiny应用程序中，可方便轻松访问，并且允许用户尝试设置不同的参数值来查看其疾病的后果。查看结果的一种方法是绘制四个群体中每个群体的人口比例：

 1library("ggplot2") 2library("tidyverse") 3ldata <- as.data.frame(out) %>% gather(key, value, -time)  4head(ldata) 5##   time key     value 6## 1  0.0   S 0.1449991 7## 2  0.2   S 0.1449991 8## 3  0.4   S 0.1449991 9## 4  0.6   S 0.144999010## 5  0.8   S 0.144999011## 6  1.0   S 0.144999012#plot13ggplot(data=ldata, 14       aes(x = time,15           y = value,16           group = key,17           col = key18           )) + 19      ylab("Proportion of full population") + xlab("Time (days)") +20      geom_line(size = 2) + 21      scale_colour_manual(values = c("red", "green4", "black", "blue")) +22      scale_y_continuous(labels = scales::percent, limits = c(0, 1))

该图显示了传染病模型的结果以及随时间各人群所占的比例。有趣的是观察到的易感性比例下降的程度（以及多快）。蓝线显示接种疫苗/先前免疫个体的实际比例，在该模型中它一直保持不变。当疾病暴发时，恢复组中的新个体将被进入到免疫组中。还可以计算显示爆发影响的简单统计数据。

1# Proportion of full population that got the 2# disease by end of time frame3ldata %>% filter(time == max(time), key=="R") %>% 4          mutate(prop = round(100 * value, 2))5## Warning: package 'bindrcpp' was built under R version 3.5.26##   time key     value  prop7## 1  200   R 0.1137451 11.37

在这里，我们可得出当疾病暴发时，约11.37％的人口患有此病。我们还可以计算患有疾病的易感人群的比例。

1# Proportion of susceptible that will get 2# the disease by end of time frame3as.data.frame(out) %>% filter(row_number() == n()) %>% 4          mutate(res = round(100*(R + I) / (S + I + R), 2)) %>% pull(res)5## [1] 81.84

如何预防疾病的流行

R0表示受感染者将感染的人数。如果这个数字小于1，那么爆发将自行消失，如果R0> 1则爆发流行病。群体免疫力是接种/免疫人数所占的百分比，与R0密切相关。我们可以直接计算这些数字：

1# Proportion of population that need to be 2# immune for herd immunity3round(100 * (1 - 1 / (connum)), 2)4## [1] 93.33

因此，当R0 = 15时，从第1天开始需要93.33％的免疫力。第1天的有效R0是根据已经免疫的人数按比例缩小的疾病的R0。尽管这里的基本生殖数是R0 = 15，但由于免疫和接种疫苗的个体数量，感染者感染的个体的有效数量仅为2.18。但是，由于这个数字大于1，它仍然会变成爆发。

在评估疫苗接种的影响时，需要考虑与接种疫苗的成本和风险相比的疾病发生的频率，风险和成本。例如，假设1000名感染麻疹的人中大约有1-2人会死亡，10000名接种疫苗的群体中如有1人面临大脑水肿的风险。如果某一地区570万人中只有1％的人感染麻疹，那么对应的麻疹病人就有57000人，其中%2的死亡率估计，114人可能因麻疹而死亡，但是按照80%接种率，接种疫苗后会有456人将面临大脑水肿的风险。

应该强调的是，我们已经使用常微分方程进行所有计算。因此，我们没有考虑在实际情况中可能发生的随机波动。因此，我们在这里看到的数字应该被视为平均效应，而实际效果可能比我们在这里看到的更温和或更差。

传染病模型Shiny应用

源代码获取：https://github.com/ekstroem/Shiny-Vaccine

文章来源： Buildinga Shiny app to show the impact of vaccines

http://sandsynligvis.dk/2019/03/06/building-a-shiny-app-to-show-the-impact-of-vaccines/

往期精彩：

• 诹图系列(1): 简单条形图

• 字节跳动(今日头条)，为何战斗力如此凶猛？

• shinydashboard与shiny_史上最全（二）

• R语言中文社区2018年终文章整理（作者篇）

• R语言中文社区2018年终文章整理（类型篇）