武汉加油，中国加油--全国新型肺炎计算记录(R0,SIR)

武汉加油，中国加油基本再生数**基本再生数（basic reproduction number，R0）**是指在一个全是易感染态个体构成的群体中，一个感染态的个体在恢复之前平均能感染的人数。在流行病学中，R0>1 表示疾病将爆发，<1 则表示疾病走向消亡，故 R0 是判断流行病是否爆发的重要条件之一。但在真实疫情的传播过程中，由于政府干预政策实施、个体行为改变（戴口罩、减少出行等）..

雪球干死黄旭东

10212人浏览 · 2020-02-08 12:25:07

雪球干死黄旭东 · 2020-02-08 12:25:07 发布

武汉加油，中国加油

基本再生数

**基本再生数（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），指被隔离，或因病愈而具有免疫力的人。
来自百度百科

数据来源

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之间.