武汉加油,中国加油
基本再生数
**基本再生数(basic reproduction number,R0)**是指在一个全是易感染态个体构成的群体中,一个感染态的个体在恢复之前平均能感染的人数。在流行病学中,R0>1 表示疾病将爆发,<1 则表示疾病走向消亡,故 R0 是判断流行病是否爆发的重要条件之一。但在真实疫情的传播过程中,由于政府干预政策实施、个体行为改变(戴口罩、减少出行等)、易感人群数量减少(因患病人数增多或使用疫苗等)等外在因素影响,不再满足基本再生数定义的理想模型条件。为刻画流行病随时间的演化过程,我们将传播过程中某一时刻下(t)平均一个感染态个体感染的人数定义为有效再生数,记为 Rt。在实际疫情的控制过程中,当 Rt<1 时,即平均一个感染态个体感染的人数小于 1 时,认为疾病已被控制同时疾病也将走向消亡。
SIR模型
传染病模型有着悠久的历史,一般认为始于1760年Daniel Bernoulli在他的一篇论文中对接种预防天花的研究。真正的确定性传染病数学模型研究的前进步伐早在20世纪初就开始了,Hamer、Ross等人在建立传染病数学模型的研究中做出了大量的工作,直到1927年Kermack与McKendrick在研究流行于伦敦的黑死病时提出了的SIR仓室模型,并于1932年继而建立了SIS模型,在对这些模型的研究基础上提出了传染病动力学中的阈值理论。Kermack与McKendrick的SIR模型是传染病模型中最经典、最基本的模型,为传染病动力学的研究做出了奠基性的贡献。
模型中把传染病流行范围内的人群分成三类:S类,易感者(Susceptible),指未得病者,但缺乏免疫能力,与感病者接触后容易受到感染;I类,感病者(Infective),指染上传染病的人,它可以传播给S类成员;R类,移出者(Removal),指被隔离,或因病愈而具有免疫力的人。
来自百度百科
相关论文
基本再生数
武汉新型冠状病毒感染肺炎基本再生数的初步预测
论文地址:http://www.cjebm.com/article/10.7507/1672-2531.202001118
对源自中国武汉的2019-nCoV爆发的潜在国内和国际传播的预测和预测:一项模型研究
论文地址:ttps://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30260-9/fulltext
武汉市2019例新型冠状病毒性肺炎99例流行病学和临床特征
论文地址:https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30211-7/fulltext
与2019年新型冠状病毒相关的家族性肺炎,表明人与人之间的传播:家庭聚类研究
论文地址:https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30154-9/fulltext
中国武汉市2019年新型冠状病毒感染患者的临床特征
论文地址:https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30183-5/fulltext
SIR模型
《SIR模型在成人麻疹爆发及其疫情控制评价中的应用》山东大学学报(医学版)2016年 第54卷 第9期
《H1N1 甲型流感全球航空传播与早期预警研究》《中国科学 》杂志社 2010 年 第55卷 第12期
数据来源
https://voice.baidu.com/act/newpneumonia/newpneumonia/?from=osari_pc_3
以及每日深圳公布的详细感染病例
代码部分
# from login.mysql_login import conn as CONN
import pandas as pd,numpy as np
from scipy.integrate import odeint
import math,datetime
import re
def getCureTime(cure_time:str):
cure_time = cure_time.split('。')[-2]
cure_time = cure_time.replace('已','')
cure_time = cure_time.replace('痊愈出院','')
cure_time = cure_time.replace('治愈出院','')
cure_time = re.sub('(.+?)','',cure_time)
if cure_time == '目前情况稳定':return None
if cure_time == '23日':
cure_time = '1月' + cure_time
cure_time = '2020年' + cure_time
cure_time = datetime.datetime.strptime(cure_time, "%Y年%m月%d日")
return cure_time
if __name__ == '__main__':
# sp_ncov_status_all = pd.read_sql('''
# select occur_period_f,
# admin_area_name,
# confirm as '确诊病例(人)',
# cure as '治愈人数(人)',
# death as '死亡人数(人)',
# suspect as '疑似病例(人)'
# from qhdata_support_spider.sp_ncov_status
# WHERE admin_area_code = 0086
# ORDER BY occur_period_f
# ''',CONN)
#
# sp_ncov_status_sz = pd.read_sql('''
# select occur_period_f,
# admin_area_name,
# confirm as '确诊病例(人)',
# cure as '治愈人数(人)',
# death as '死亡人数(人)',
# suspect as '疑似病例(人)'
# from qhdata_support_spider.sp_ncov_status
# WHERE admin_area_code = 440300000000
# ''', CONN)
sp_ncov_status_all = pd.read_excel('source/sp_ncov_status_all.xlsx')
sp_ncov_status_sz = pd.read_excel('source/sp_ncov_status_sz.xlsx')
detail_case = pd.read_excel('source/深圳市冠状病毒病例个案数据(汇总)-20200211.xlsx')
# ##############################数据预处理##############################################
for i in ['确诊病例(人)','治愈人数(人)','死亡人数(人)']:
# sp_ncov_status_all['当日' + i] = sp_ncov_status_all[i] - sp_ncov_status_all[i].shift()
sp_ncov_status_sz['当日' + i] = sp_ncov_status_sz[i] - sp_ncov_status_sz[i].shift()
detail_case['疑似到确诊时间间隔'] = detail_case[['发病时间','确诊时间']].\
apply(lambda x:None if x[0] == '-' else
None if x[1] == '-' else
(x[1] - x[0]).days,
axis=1)
detail_case['疑似到确诊时间间隔'] = detail_case['疑似到确诊时间间隔'].apply(lambda x:None if x < 0 else x)
detail_case['年龄'] = detail_case['年龄'].apply(lambda x:x.replace('。','') if type(x) == str else x)
detail_case['年龄'] = detail_case['年龄'].astype(int)
detail_case['治愈时间'] = detail_case[['目前是否治愈','完整病例介绍原文']]\
.apply(lambda x:getCureTime(x[1]) if x[0] == '是' else None,
axis=1)
detail_case['恢复天数'] = detail_case[['确诊时间','治愈时间']].\
apply(lambda x:None if x[0] == '-' else
None if x[1] == '-' else
(x[1] - x[0]).days,
axis=1)
# 以第一个不明原因肺炎发现者的报道时间2019年12月8日为感染的第一天
first_date = datetime.datetime.strptime('2019-12-08', "%Y-%m-%d")
sp_ncov_status_all['t'] = sp_ncov_status_all['occur_period_f'].apply(lambda x:(datetime.datetime.strptime(x, "%Y-%m-%d") - first_date).days + 1)
# #####################################################################################
# ##############################数据统计################################################
case_num = detail_case.shape[0]
# 性别比例
sex_info = pd.pivot_table(detail_case,
index=['性别'],
values=['病例序号'],
aggfunc={'病例序号': [len,lambda x:len(x)/case_num]}).reset_index()
sex_info.columns = ['性别','占比','人数']
# 年龄分层
temp_df = detail_case[['年龄','病例序号']]
temp_df['年龄分组'] = temp_df.年龄.apply(lambda x:'1.青少幼年(20以下)' if x < 20 else
'2.青壮年(20-50)' if ((x >= 20) and (x < 50)) else
'3.中年(50-70)' if ((x >= 50) and (x < 70)) else
'4.老年(70以上)'
)
age_info =temp_df[['年龄分组','病例序号']].\
groupby('年龄分组').\
agg({'病例序号':[len,lambda x:len(x)/case_num]}).reset_index()
age_info.columns = ['年龄分组','人数','占比']
# 疑似到确诊时间间隔
diagnostic_interval = pd.DataFrame(detail_case.describe()['疑似到确诊时间间隔']).T
diagnostic_interval_1 = pd.pivot_table(detail_case,
index=['疑似到确诊时间间隔'],
values=['病例序号'],
aggfunc={'病例序号': [len,lambda x:len(x)/case_num]}).reset_index()
diagnostic_interval_1.columns = ['疑似到确诊时间间隔','占比','人数']
# 恢复天数
restore_day = pd.DataFrame(detail_case.describe()['恢复天数']).T
restore_day_1 = pd.pivot_table(detail_case,
index=['恢复天数'],
values=['病例序号'],
aggfunc={'病例序号': [len,lambda x:len(x)/case_num]}).reset_index()
restore_day_1.columns = ['恢复天数','占比','人数']
# #####################################################################################
# ##############################R0计算逻辑##############################################
avg_day = diagnostic_interval['mean'][0].round(2) # 疑似到确诊时间间隔引用深圳市的平均值(天)
# 选择拥有新增最新数据的2.4-2.6日的确认数据与对应的2.1-2.3疑似数据(相隔4天左右),计算转化率p
p = round((3887 + 3694 + 3151) / (4562 + 5173 + 5072), 2) # 这里计算出来的p大约0.72左右
# SARS传播的p的取值在0.5-0.8之间;初期报导提到59个疑似病例有41位最终确诊,因此p的一个参考值为41/59=0.695;
# 以第一个不明原因肺炎发现者的报道时间2019年12月8日为t=1
# 2020年2月6日实际上被感染人数为Y(t),t=61
Y_t = 3151 + round(4833*p,0)
t = 61
T_g = 8.4# 生成时间(generation time)
T_i = (5+6+4+3)/4 # 潜伏时间(incubation time)
lamda = math.log(Y_t,math.e)/t
Rho = T_i/T_g # 潜伏期占生成时间的比例Rho
# 论文里面说平均潜伏期5.1天,不过我没有找到。钟南山说3天。
R_0 = round(1 + lamda*T_g + Rho*(1 - Rho) * math.pow((lamda*T_g),2),2)
# #####################################################################################
# ##############################R0数据计算##############################################
def basicReproductionNumber(Y_t,t,T_i,T_g = 8.4):
lamda = math.log(Y_t,math.e)/t
Rho = T_i/T_g # 潜伏期占生成时间的比例Rho
R_0 = 1 + lamda * T_g + Rho * (1 - Rho) * math.pow((lamda * T_g), 2)
return R_0
sp_ncov_status_all_temp = sp_ncov_status_all.dropna()
sp_ncov_status_all_temp['Y_t'] = sp_ncov_status_all_temp['新增确诊(人)'] + sp_ncov_status_all_temp['新增疑似(人)'] * p
sp_ncov_status_all_temp['Y_t'] = sp_ncov_status_all_temp['Y_t'].round(0)
sp_ncov_status_all_temp['R0'] = sp_ncov_status_all_temp[['Y_t','t']].\
apply(lambda x:basicReproductionNumber(x[0],x[1],T_i,8.4),axis=1)
# #####################################################################################
# ##############################SIR计算逻辑################################################
# https://towardsdatascience.com/modelling-the-coronavirus-epidemic-spreading-in-a-city-with-python-babd14d82fa2
# SIR模型中的总人口数, N.
N = 11081000 # 这里采用2019年武汉市常住人口数(百度的)
# 感染者I与移除者R
I0, R0 = 44, 0 # I0选择最初的武汉华南海鲜市场的44个感染者作为最初的感染者
# 易感人群S
S0 = N - I0 - R0
# 恢复率(移除率)gamma
# 这里是选取深圳感染案例的平均恢复天数,实际上不同地区的医疗水平不同会导致不同的恢复天数。
gamma = 1. / restore_day['mean'][0]
# 接触率beta
# 这里的R_0是全国的R_0,但是实际上应该将R_0分成湖北省内与全国除湖北省的两种情况。
# 现在省内R_0明显比全国的要高的多,而省外明显要比全国的低的多。在此我们额外添加数个R_0来观测[0.9,2,3,4]
# R_0 = 4
beta = R_0 * gamma # 基本再生数 R0, 等于接触率除以移出率
# 最长时间(天)
day = 250
t = np.linspace(0, day, day)
def deriv(y, t, N, beta, gamma):
S, I, R = y
dSdt = -beta * S * I / N
dIdt = beta * S * I / N - gamma * I
dRdt = gamma * I
return dSdt, dIdt, dRdt
y0 = S0, I0, R0
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T
# #####################################################################################
# ##############################画图###################################################
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
fig = plt.figure(figsize=[16, 9]) # 设定图片大小
temp_df = age_info[['年龄分组','人数']]
sns.barplot(x='年龄分组', y='人数',data=temp_df[['年龄分组','人数']])
plt.savefig('pic/年龄分层.png')
plt.close()
fig = plt.figure(figsize=[16, 9]) # 设定图片大小
ax1 = fig.add_subplot(111) # 添加第一副图
ax2 = ax1.twinx() # 共享x轴,这句很关键!
sns.barplot(x='性别', y='人数', data=sex_info[['性别','人数']], ax=ax1, alpha=0.8) # 画柱状图,注意ax是选择画布
sns.pointplot(x='性别', y='占比', data=sex_info[['性别','占比']], ax=ax2) # 画在ax2画布上
plt.savefig('pic/性别比例.png')
# plt.show()
plt.close(fig)
fig = plt.figure(figsize=[16, 9]) # 设定图片大小
temp_df = diagnostic_interval_1[['疑似到确诊时间间隔','人数']]
sns.barplot(x='疑似到确诊时间间隔', y='人数',data=temp_df[['疑似到确诊时间间隔','人数']])
plt.savefig('pic/疑似到确诊时间间隔.png')
plt.close()
fig = plt.figure(figsize=[16, 9]) # 设定图片大小
temp_df = restore_day_1[['恢复天数','人数']]
sns.barplot(x='恢复天数', y='人数',data=temp_df[['恢复天数','人数']])
plt.savefig('pic/恢复天数.png')
plt.close()
fig = plt.figure(figsize=[16, 9]) # 设定图片大小
temp_df = sp_ncov_status_all_temp[['occur_period_f', 'R0']]
sns.pointplot(x='occur_period_f', y='R0',data=temp_df)
plt.savefig('pic/R0变化趋势.png')
plt.close()
# SIR模型
fig = plt.figure(facecolor='w',figsize=[16, 9],dpi=144)
plt.style.available # style样式 https://neusncp.com/user/blog?id=203
plt.style.use('seaborn-darkgrid')
plt.suptitle('R_0: %s' % R_0, fontsize=16, fontweight='bold')
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
ax = fig.add_subplot(111, axisbelow=True)
ax.plot(t, S / 10000, 'b', alpha=0.5, lw=2, label='易感人群')
ax.plot(t, I / 10000, 'r', alpha=0.5, lw=2, label='感染者')
ax.plot(t, R / 10000, 'g', alpha=0.5, lw=2, label='移除者')
ax.set_xlabel('时间(天)')
ax.set_ylabel('人数(万人)')
ax.xaxis.set_major_locator(plt.MultipleLocator(10))
ax.set_xlim(0, day)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(False)
plt.savefig(f'pic/SIR变化趋势(R_0={R_0}).png')
plt.close()
# #####################################################################################
部分说明
1.计算公式的来由:
中国循证医学杂志的《武汉新型冠状病毒感染肺炎基本再生数的初步预测》
柳叶刀的《Nowcasting and forecasting the potential domestic and
international spread of the 2019-nCoV outbreak originating
in Wuhan, China: a modelling study》
2.Tg的来由:
柳叶刀的《A familial cluster of pneumonia associated with the 2019
novel coronavirus indicating person-to-person transmission:
a study of a family cluster》
红框是没有去武汉的深圳一家人中的一人,从中可以看出,该人在其他家庭成员回来后大约经过8.4天后发病。出现初步症状是在8号左右。
3.Tg的来由:
柳叶刀的《A familial cluster of pneumonia associated with the 2019
novel coronavirus indicating person-to-person transmission:
a study of a family cluster》
不过在论文《武汉新型冠状病毒感染肺炎基本再生数的初步预测》(http://www.cjebm.com/article/10.7507/1672-2531.202001118)中写的平均潜伏时间是5.1天,具体我在源论文中没有找到。我只找到这个表格有个类似的数据。
结果展示
年龄分组
性别分组
疑似到确诊时间间隔分布
恢复天数
全国R0变化趋势
简易SIR变化趋势
这里的t是以第一个不明原因肺炎发现者的报道时间2019年12月8日为感染的第一天。
SIR模型中的移除者可以使治愈,可以是隔离,也可以是死亡,请注意。
从由于R_0采用的是全国数据来计算的,但是实际上应该将R_0分成湖北省与全国除湖北省两种情况来看待。从现在看全国(除湖北)根本没有爆发大规模感染,所以我假设全国(除湖北)的R_0小于1。而湖北到今日(2月14日,与2019年12月8日间隔69天)其确诊、疑似病患数量仍然上升,可以看出其R_0比全国R_0(2.58)要高,所以这里使用3,4来测试。
从图与实际结合推测,全国(除湖北)的SIR比较接近于R_0=0.9情况,而湖北则从现在仍然是爆发期来看,比较接近R_0=3-4之间.