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

武汉加油，中国加油

基本再生数

**基本再生数（basic reproduction number，R0）**是指在一个全是易感染态个体构成的群体中，一个感染态的个体在恢复之前平均能感染的人数。在流行病学中，R0>1 表示疾病将爆发，<1 则表示疾病走向消亡，故 R0 是判断流行病是否爆发的重要条件之一。但在真实疫情的传播过程中，由于*干预政策实施、个体行为改变（戴口罩、减少出行等）、易感人群数量减少（因患病人数增多或使用疫苗等）等外在因素影响，不再满足基本再生数定义的理想模型条件。为刻画流行病随时间的演化过程，我们将传播过程中某一时刻下（t）平均一个感染态个体感染的人数定义为有效再生数，记为 Rt。在实际疫情的控制过程中，当 Rt<1 时，即平均一个感染态个体感染的人数小于 1 时，认为疾病已被控制同时疾病也将走向消亡。

相关论文

武汉新型冠状病毒感染肺炎基本再生数的初步预测
论文地址: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

数据来源

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
import math,datetime




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/深圳市冠状病毒病例个案数据（汇总）-20200206.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['年龄'].astype(int)

    # 以第一个不明原因肺炎发现者的报道时间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 = ['疑似到确诊时间间隔','占比','人数']
    # #####################################################################################
    # ##############################计算逻辑################################################
    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)
    # 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）
    lamda = math.log(Y_t,math.e)/t
    Rho = (T_g - avg_day)/T_g # 潜伏期占生成时间的比例Rho
    R_0 = 1 + lamda*T_g + Rho*(1 - Rho) * math.pow((lamda*T_g),2)
    # #####################################################################################
    # ##############################数据计算################################################
    def basicReproductionNumber(Y_t,t,avg_day,T_g  = 8.4):
        lamda = math.log(Y_t,math.e)/t
        Rho = (T_g - avg_day) / 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],avg_day,8.4),axis=1)
    # #####################################################################################
    # ##############################画图###################################################
    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 = 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()

# #####################################################################################

部分说明

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号左右。

结果展示

年龄分组

性别分组

疑似到确诊时间间隔分布

R0变化趋势