一、引言
随着时间的推移，新冠病毒的二次感染人群越来越多，这意味着在春节过后暂时消失的病毒又有可能卷土重来。目前，我国主流的病毒株为XBB，中疾控中心的数据显示，自五月上旬到五月中旬，本土病例中XBB的占比已经从83.6%增加到91.9%。
近期医院发热门诊人数增多的情况也表明第二波疫情可能已经到来。钟南山院士在5月22日的发言中指出，第二波疫情高峰很可能会在6月底到来。现在已经进入6月，如果院士的预测没有出错，那么二阳的人数将在未来的日子里呈现大幅度增加的趋势。
图中可以看出，「终南山院士使用的预测模型是seirs，今天我们来研究和学习一下seirs是什么？如何使用？」
二、seirs简介
「SEIRS模型」是一种流行病学模型，用于「建模呼吸道传染病的传播」。在SEIRS模型中，人口根据其在传染病生命周期中的感染状态被划分为四个类别：「易感（susceptible）、暴露（exposed）、感染（infectious）和康复（recovered）」。SEIRS模型通过考虑个体的人群动态和传染病传播的基本机制，对传染病疫情的演化进行建模。下面是该模型各个状态的定义及其数学表示：
- 「易感状态（S）」：未感染病毒的人群，可以被感染。
- 「暴露状态（E）」：感染病毒但还未出现症状的人群。在这个状态下，人群处于潜伏期，也称为隐藏状态。
- 「感染状态（I）」：患病且具有传染性的人群。
- 「康复状态（R）」：已从疾病中康复的人群，对后续的疾病传播没有贡献。
在SEIRS模型中，每个类别的人群都可以以以下方式变化：
- 「易感状态」：易感人群可以转变为暴露状态，通过感染的机会是易感人群所接触到已确诊感染的人群。
- 「暴露状态」：暴露人群可以转变为感染状态。根据COVID-19的传播方式，一般潜伏期在9天左右，在这期间，暴露人群会变得有传染性，这也是模型中一个重要的变量。
- 「感染状态」：感染人群可以转变为康复状态，也可以由于治愈失败而转变为死亡状态。它们所感染的人群一般是其他易感人群，传染性与病毒传染的强度有很大关系。 有一部分从疾病中恢复的人群变成康复状态，对此我将其称为R群体。这个状态表示患病者从病毒中恢复，具备免疫力。
- 「康复状态」：康复状态的人群一般来说可以认为会对病毒具有抵御力，因此在后续的传染爆发阶段中可以不作为传染源。
SEIRS模型可以用以下一系列微分方程来表示：
dS/dt= -beta * S * I / N dE/dt= beta * S * I / N - sigma * E dI/dt= sigma * E - gamma * I dR/dt= gamma * I
其中，「N」是总人口，即 「N=S+E+I+R」，即易感、隐藏、感染、康复人群的总和，而 「beta」表示感染概率。潜伏期的时间长度与 「sigma」 相关，「gamma表示康复率」。这些参数根据特定的传染病和特定的人群情况而有所不同。您可以根据现实中的数据和情况来调整这些参数，以便更准确地模拟传染病的传播，并预测其可能的发展趋势。
三、数据来源
根据目前的官方数据（参考国家卫生健康委员会的官方网站http://www.nhc.gov.cn/），新冠病毒在中国的传播率（传染率）大约为0.25，即每名感染者平均可以传播给0.25名未感染者，而治愈率（包括确诊和无症状感染者）达到了约99.5%，表明大部分感染者都能够康复和治愈。此外，新冠病毒的平均潜伏期在中国约为14天，即感染病毒后，患者在14天内可能仍然没有明显症状，但仍可以对他人进行传染。值得注意的是，这些数据仅供参考，实际情况会根据时间和地点的变化而有所不同。即：
install.packages("deSolve") library(deSolve) N <- 1400000000 # 中国总人口数量 beta <- 0.13 # 病毒传播速率 gamma <- 0.995 # 治愈率（包括确诊与无症状感染者） sigma <- 1/14 # 平均潜伏期 I <- 41699 #初始感染人数
四、处理方法
library(deSolve) # 定义 SEIRS 模型 seirs <- function(time, state, parms) { with(as.list(c(state, parms)), { # 模型方程 dS <- -beta*S*I/(S+E+I+R) dE <- beta*S*I/(S+E+I+R) - rho*E dI <- rho*E - gamma*I dR <- gamma*I # 返回微分方程组 return(list(c(dS, dE, dI, dR))) }) } # 设置模型参数 parms <- c(S = 1.4e9, E = 1, I = 41699, R = 0, rho = 1/14, gamma = 0.13, beta = 0.25) # 设置初始状态 init <- c(S = 1.4e9, E = 1, I = 41699, R = 0) # 设置起始时间为今天 start_time <- as.Date("2022-12-23") # 模拟未来 N 天的趋势 N <- 365 times <- seq(start_time, by = "day", length.out = N) times <- times <- as.numeric(times) out <- as.data.frame(ode(y = init, times = times, func = seirs, parms = parms)) # 可视化结果 library(ggplot2) df <- data.frame(time=as.Date(out$time, origin = "1970-01-01") , S=out[,2], E=out[,3], I=out[,4], R=out[,5]) # 美化图表及设置线条属性 theme_set(theme_minimal(base_size = 14)) # 设置主题 line_size <- 1.5 # 设置线的宽度 line_type <- c('solid', 'dotted', 'dashed', 'dotdash') # 设置不同种类线条的样式 # 绘制图表 ggplot(df, aes(x = time)) + geom_line(aes(y = S, color = "Susceptible"), size = line_size, linetype = line_type[1]) + geom_line(aes(y = E, color = "Exposed"), size = line_size, linetype = line_type[2]) + geom_line(aes(y = I, color = "Infected"), size = line_size, linetype = line_type[3]) + geom_line(aes(y = R, color = "Recovered"), size = line_size, linetype = line_type[4]) + scale_color_manual(values = c("blue", "orange", "red", "green")) + # 自定义颜色 labs(title = "SEIR Model for COVID-19 in China", x = "Time (days)", y = "Number of individuals", color = "Cases") + theme(legend.position = "bottom", legend.title = element_blank(), axis.title = element_text(size = 16), axis.text = element_text(size = 12), plot.title = element_text(size = 18, face = "bold", hjust = 0.5), panel.border = element_rect(colour = "black", fill = NA), panel.grid.major = element_line(colour = "gray", linetype = "dotted", size = 0.5), panel.grid.minor = element_line(colour = "gray", linetype = "dotted", size = 0.3)) # 输出峰值时间和数量 peak_time <- df$time[which.max(df$I)] peak_count <- max(df$I) cat(paste("Peak time: ", peak_time, " days")) cat("\n") cat(paste("Peak count: ", peak_count, " individuals"))
结果展示：
Peak time: 2023-09-22 Peak count: 69027262.4873488
五、结果解释
我们的算出的结果疫情高发期是2023-09-22和峰值数量是69027262。和「终南山院士」算出的时间出入了2个月左右，然后峰值数量差距400多万，还是差距很大的。我获取的数据不够准确全面，然后没有加入过去的案例辅助预测，可能初始参数设置有差别和数据质量的原因导致了巨大的偏差。
「在使用 SEIRS 模型预测传染病的传播时，正确地设置初始参数十分重要。例如，初始感染人数、初始接触人数等参数可能会对结果产生很大的影响。因此，需要对模型进行敏感性分析，确定最优的初始参数。」
