如何简单构建新冠肺炎的预测模型?——附R、python、matlab代码

专题合集更多教程

作者:李健民

 

疫情开始的时候,我就坚守在发热门诊,下班后开始思考如何构建新冠肺炎的预测模型。查阅相关书籍和文献后,我发现这属于传染病动力学模型的范围,已经有一套成熟的理论。

 

随着疫情的进展,网上相关的文章也越来越多,但是里面的内容涉及比较多微积分的知识,对一线临床医生不太友好。最近正好在B站看见毕导的科普视频(https://www.bilibili.com/video/av85508117),我就以截图的参数为例,构建一个简单的SEIR模型,并附上R、python、matlab代码。对代码不感兴趣的同学可直接看原理和最后总结即可。

 

 

1、什么是SEIR模型

 

常见的传染病模型按照传染病类型分为 SI、SIR、SIRS、SEIR 模型等,用于研究传染病的传播速度、空间范围、传播途径、动力学机理等问题,以指导对传染病的有效地预防和控制。

 

首先介绍S、E、I、R几个重要的参数

1、S 类:易感者 (Susceptible),指未得病者,但缺乏免疫能力,与感染者接触后容易受到感染;在视频中,假设某区域的人口数为10000,那么第一天的S=N-I=9999

2、E 类:暴露者 (Exposed),指接触过感染者,但暂无能力传染给其他人的人,对潜伏期长的传染病适用;本例中第一天为0个。

3、I 类:感病者 (Infectious),指染上传染病的人,可以传播给 S 类成员,将其变为 E 类或 I 类成员;本例中第一天为1个。

4、R 类:康复者 (Recovered),指被隔离或因病愈而具有免疫力的人。如免疫期有限,R 类成员可以重新变为 S 类。本例中第一天为0个。

 

接下来看看图中的r、β、γ、α

1、r:感染患者(I)每天接触的易感者数目,本例为20

2、β:传染系数;由疾病本身的传播能力,人群的防控能力决定,本例设置为0.03

3、γ:恢复系数;一般为病程的倒数,例如流感的病程5天的话,那么它的γ就是1/5,本例设置为0.1

4、α:潜伏者的发病概率;一般为潜伏期的倒数,本例为0.1

 

2、怎么理解SEIR模型

当dt=1时

St+1-St=-γβItSt / N

Et+1-Et=γβItSt / N-ɑEt

It+1-It=ɑEt-γIt

Rt+1-Rt=γIt

 

我们看流程图和公式,dS/dt可以理解为当时间t无限接近0时,S的变化量。当t=1时,dS/dt就是每日S的改变数,其余的dE/dt、dI/dt、dR/dt同理。因为S是从一开始接近N,然后慢慢被感染成E,呈下降变化,所以第一条公式右边是带负号。

 

每天有多少S减少由每天发病人群接触人数(r)、传染系数(β)、发病人数的比例(I/N)、易感人群的比例(S/N)、总人数(N)所决定的,所以dS/dt=-r*β*S/N*I/N*N=-rβSI/N。那么每天E的改变数量dE/dt就是每天S转化为E的数目(rβSI/N)减去每天E转化为I的数目(rβSI/N-αE)。其余的dI/dt、dR/dt同理。

 

3、构建R代码

 

# 定义初始状态各人数

N = 10000

E = 0

I = 1

S = N-I

R = 0

 

# 常量赋值

r = 20

B = 0.03

a = 0.1

y = 0.1

 

#设置时间观察到第150天

T = 150

 

#用for循环的方式构建上述方程组

for (i in 1:(T-1)){

  S[i+1] = S[i] - r*B*S[i]*I[i]/N

  E[i+1] = E[i] + r*B*S[i]*I[i]/N - a*E[i]

  I[i+1] = I[i] + a*E[i] - y*I[i]

  R[i+1] = R[i] + y*I[i]

 

}

 

#生成表格并查看,表格就是每天S、E、I、R的值

result <- data.frame(S, E, I, R)

View(result)

 

#作图

X_lim <- seq(1,T,by=1)

plot(S~X_lim, pch=15, col="DarkTurquoise", main = "SEIR Model", type = "l", xlab = "时间T", ylab = "各人数", xlim = c(0,T), ylim = c(0,10000))

lines(S, col="DeepPink", lty=1) 

lines(E, col="DarkTurquoise", lty=1)

lines(I, col="RosyBrown", lty=1)

lines(R, col="green", lty=1)

legend(120,8000,c("S","E","I","R"),col=c("DeepPink","DarkTurquoise","RosyBrown"),text.col=c("DarkTurquoise","DeepPink","RosyBrown","Green"),lty=c(1,1,1,1))

 

 

查看result表格,我们可以看到在第77天的时候,感染者达到峰值

 

4、构建matlab代码

 

毕导视频的图是用matlab做的,那么这里也做一个

%--------------------------------------------------------------------------

%   初始化

%--------------------------------------------------------------------------

clear;clc;

 

%--------------------------------------------------------------------------

%   参数设置

%--------------------------------------------------------------------------

N = 10000;                                                                  

E = 0;                                                                     

I = 1;                                                                      

S = N - I;                                                                  

R = 0;                                                                     

 

r = 20;                                                                    

B = 0.03;                                                                  

a = 0.1;                                                                 

y = 0.1;                                                                   

 

T = 1:150;

for i = 1:length(T)-1

    S(i+1) = S(i) - r*B*S(i)*I(i)/N(1);

    E(i+1) = E(i) + r*B*S(i)*I(i)/N(1)-a*E(i);

    I(i+1) = I(i) + a*E(i) - y*I(i);

    R(i+1) = R(i) + y*I(i);

    

end

 

plot(T,S,T,E,T,I,T,R);grid on;

xlabel('天');ylabel('人数')

legend('易感者','潜伏者','传染者','康复者')

 

最后得出图像,可以看到R和matlab的代码基本一致,后者的作图更好看。有兴趣的同学可以试试用ggplot在R上重新作图

 

5、构建python代码

 

有少数临床的同学也会用python做数据分析,我参考网上的资料后(https://blog.csdn.net/shelgi/article/details/104108757),补充一份python的代码

 

import math

import numpy as np

import matplotlib

import matplotlib.pyplot as plt

 

plt.rcParams['font.sans-serif'] = ['KaiTi']

plt.rcParams['axes.unicode_minus'] = False

 

def calc(T):

    for i in range(0, len(T) - 1):

        S.append(S[i] - r * b * S[i] * I[i] / N )

        E.append(E[i] + r * b * S[i] * I[i] / N - a * E[i] )

        I.append(I[i] + a * E[i] - y * I[i])

        R.append(R[i] + y * I[i])

 

def plot(T,S,E,I,R):

    plt.figure()

    plt.title("SEIR-病毒传播时间曲线")

    plt.plot(T,S,color='r',label='易感者')

    plt.plot(T, E, color='k', label='潜伏者')

    plt.plot(T, I, color='b', label='传染者')

    plt.plot(T, R, color='g', label='康复者')

    plt.grid(False)

    plt.legend()

    plt.xlabel("时间(天)")

    plt.ylabel("人数")

    plt.show()

 

if __name__ == '__main__':

    # 首先还是设置一下参数,之后方便修改

    N = 10000  # 人口总数

    E = []  # 潜伏携带者

    E.append(0)

 

    I = []  # 传染者

    I.append(1)

 

    S = []  # 易感者

    S.append(N - I[0])

 

    R = []  # 康复者

    R.append(0)

 

    r = 20  # 传染者接触人数

    b = 0.03  # 传染者传染概率

    a = 0.1  # 潜伏者患病概率

    y = 0.1  # 康复概率

 

    T = [i for i in range(0, 160)]  # 时间

    calc(T)

    plot(T,S,E,I,R)

 

最后作图如下

 

6、总结与不足

 

本例只作为简单的SEIR入门介绍,现实情况肯定要复杂很多,譬如毕导视频后面基于潜伏期患者也能传播这一事实,对SEIR模型做了小小的修改。

 

另外,现实中还有很多问题需要考虑,例如:如果感染者分为隔离和非隔离群体;如果潜伏期患者分为疑似感染者和疑似非感染者;如何通过已有的数据反推参数等等。目前的论文有不少就是基于这个模型而写的,例如这篇《A mathematical model for simulating the transmission of Wuhannovel Coronavirus》就在SEIR的基础上增加了蝙蝠、宿主概念。

 

 

毕导在视频中说他们大学三年级就会做这个模型,身为临床医生的我却在研究生后才学会,实在汗颜。这说明了我们临床专业在数据分析上,无论视野和能力都非常欠缺。临床上许多问题都需要转化为具体的数学建模问题,如果临床医生掌握好数学建模的能力,将会有更广阔的思维去研究临床问题。

 

扫码关注“医咖会”公众号,及时获取最新统计教程!

描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
我要提问
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
    提交问题