python分析新冠肺炎_如何简单构建新冠肺炎的预测模型？——附R、python、matlab代码...

最新推荐文章于 2023-04-18 16:18:20 发布

weixin_39790168

最新推荐文章于 2023-04-18 16:18:20 发布

阅读量3k

点赞数 4

原标题：如何简单构建新冠肺炎的预测模型？——附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

View(result)

#作图

X_lim

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;

fori = 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的基础上增加了蝙蝠、宿主概念。

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

责任编辑：