03_R0_estimation-WLD-02b.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# coding: utf-8\n",
    "\n",
    "\"\"\"\n",
    "title: R0 estimation in the World\n",
    "date: 2020/03/22\n",
    "author: okimebarun\n",
    "url: https://github.com/okimebarun/01_COVID19_analysis\n",
    "url: https://qiita.com/oki_mebarun\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import datetime\n",
    "import locale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "font = {'family' : 'IPAexGothic'}\n",
    "plt.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##########################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def makeCalcFrame(days):\n",
    "    t_1 = pd.Timestamp(2020,1,24) # 計算開始日\n",
    "    td = pd.Timedelta('1 days')\n",
    "    #\n",
    "    npd = [[t_1 + td * i, 0, 0, 0 ] for i in range(0,days)]\n",
    "    df1 = pd.DataFrame(npd)\n",
    "    df1.columns = ['date', 'Ppre','Pat', 'R0']\n",
    "    #\n",
    "    return df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mergeCalcFrame(df1, df2):\n",
    "    return pd.merge(df1, df2, on='date', how='left').fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calcR0(df, keys):\n",
    "    lp = keys['lp']\n",
    "    ip = keys['ip']\n",
    "    nrow = len(df)\n",
    "    getP = lambda s: df.loc[s, 'P'] if s < nrow else np.NaN\n",
    "    for t in range(nrow):\n",
    "        df.loc[t, 'Ppre'] = sum([ getP(s) for s in range(t+1, t + ip + 1)])\n",
    "        df.loc[t, 'Pat' ] = getP(t + lp + ip)\n",
    "        if df.loc[t, 'Ppre'] > 0:\n",
    "            df.loc[t, 'R0'  ] = ip * df.loc[t, 'Pat'] / df.loc[t, 'Ppre']\n",
    "        else:\n",
    "            df.loc[t, 'R0'  ] = np.NaN\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def showResult(df, title):\n",
    "    # R0=1\n",
    "    ptgt = pd.DataFrame([[df.iloc[0,0],1],[df.iloc[len(df)-1,0],1]])\n",
    "    ptgt.columns = ['date','target']\n",
    "    # show R0\n",
    "    plt.rcParams[\"font.size\"] = 12\n",
    "    ax = df.plot(title=title,x='date',y='R0', figsize=(10,7))\n",
    "    ptgt.plot(x='date',y='target',style='r--',ax=ax)\n",
    "    ax.grid(True)\n",
    "    ax.set_ylim(0,)\n",
    "    plt.show()\n",
    "    fig = ax.get_figure()\n",
    "    fig.savefig(\"R0_{}.png\".format(title))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##########################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def readCsvOfWorldArea(area : None):\n",
    "    # 下記URLよりダウンロード\n",
    "    # https://hackmd.io/@covid19-kenmo/dataset/https%3A%2F%2Fhackmd.io%2F%40covid19-kenmo%2Fdataset\n",
    "    fcsv = u'World-COVID-19.csv'\n",
    "    df = pd.read_csv(fcsv, header=0, encoding='sjis', parse_dates=[u'日付'])\n",
    "    # 日付, 対象国を抽出\n",
    "    if area is not None:\n",
    "        df1 = df.loc[:,[u'日付',area]]\n",
    "    else:\n",
    "        df1 = df.loc[:,[u'日付',u'世界全体の感染者']]        \n",
    "    df1.columns = ['date','Psum']\n",
    "    ## 累積⇒日次に変換\n",
    "    df2 = df1.copy()\n",
    "    df2.columns = ['date','P']\n",
    "    df2.iloc[0,1] = 0\n",
    "    ## 文字列⇒数値\n",
    "    getFloat = lambda e: float('{}'.format(e).replace(',',''))\n",
    "    ## 差分計算\n",
    "    for i in range(1,len(df1)):\n",
    "        df2.iloc[i, 1] = getFloat(df1.iloc[i, 1]) - getFloat(df1.iloc[i-1, 1] )\n",
    "    ##\n",
    "    return df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def R0inWorldArea(area, bShow:True):\n",
    "    keys = {'lp':5, 'ip':8 }\n",
    "    df1 = makeCalcFrame(60) # 60 days\n",
    "    df2 = readCsvOfWorldArea(area)\n",
    "    df = mergeCalcFrame(df1, df2)\n",
    "    df = calcR0(df, keys)\n",
    "    if bShow:\n",
    "        showResult(df, u'COVID-19 R0 ({})'.format(area))\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getWorldAreaList():\n",
    "    fcsv = u'World-COVID-19.csv'\n",
    "    df = pd.read_csv(fcsv, header=0, encoding='sjis', parse_dates=[u'日付'])\n",
    "    arealist = df.columns[1:]\n",
    "    return arealist\n",
    "\n",
    "arealist = getWorldAreaList()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def showResult2(ax, df, area):\n",
    "    # show R0\n",
    "    plt.rcParams[\"font.size\"] = 12\n",
    "    df1 = df.rename(columns={'R0':area})\n",
    "    df1.plot(x='date',y=area, ax=ax)\n",
    "\n",
    "def showResult3(dflist, title):\n",
    "    # R0=1\n",
    "    dfs = dflist[0][0]\n",
    "    ptgt = pd.DataFrame([[dfs.iloc[0,0],1],[dfs.iloc[len(dfs)-1,0],1]])\n",
    "    ptgt.columns = ['date','target']\n",
    "    ax = ptgt.plot(title='COVID-19 R0', x='date',y='target',style='r--', figsize=(10,8))\n",
    "    ax.set_yscale(\"symlog\", linthreshy=1)\n",
    "    #\n",
    "    for df, label in dflist:\n",
    "        showResult2(ax, df, label)\n",
    "    #\n",
    "    ax.grid(True)\n",
    "    ax.set_ylim(0,)\n",
    "    plt.show()\n",
    "    fig = ax.get_figure()\n",
    "    fig.savefig(\"R0_{}.png\".format(title))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#　\n",
    "title = '爆発的感染が観測された地域'\n",
    "arealist = ['中国本土の感染者数', 'イタリア', '米国', 'スペイン', 'イラン','韓国']\n",
    "dflist = [ [R0inWorldArea(area, False), area] for area in arealist]\n",
    "showResult3(dflist, title)\n",
    "#\n",
    "title = 'ヨーロッパ'\n",
    "arealist = ['イタリア',  'スペイン', 'ドイツ', 'フランス',  'スイス', '英国','オーストリア',\n",
    "            'ベルギー','ノルウェー','スウェーデン','デンマーク','ポルトガル']\n",
    "dflist = [ [R0inWorldArea(area, False), area] for area in arealist]\n",
    "showResult3(dflist, title)\n",
    "#\n",
    "title = 'アジア周辺で比較的感染が抑制されている地域'\n",
    "arealist = [   '台湾','日本','香港','シンガポール']\n",
    "dflist = [ [R0inWorldArea(area, False), area] for area in arealist]\n",
    "showResult3(dflist, title)\n",
    "#\n",
    "title = '今後感染拡大が懸念される地域'\n",
    "arealist = [  '米国','豪州','マレーシア', 'ブラジル','トルコ', 'イスラエル','インドネシア']\n",
    "dflist = [ [R0inWorldArea(area, False), area] for area in arealist]\n",
    "showResult3(dflist, title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#\n",
    "title = 'ヨーロッパ+推定'\n",
    "arealist = ['イタリア',  'スペイン', 'ドイツ', 'フランス',  'スイス', '英国','オーストリア',\n",
    "            'ベルギー','ノルウェー','スウェーデン','デンマーク','ポルトガル']\n",
    "dflist = [ [R0inWorldArea(area, False), area] for area in arealist]\n",
    "#\n",
    "# R0=1\n",
    "dfs = dflist[0][0]\n",
    "ptgt = pd.DataFrame([[dfs.iloc[0,0],1],[dfs.iloc[len(dfs)-1,0],1]])\n",
    "ptgt.columns = ['date','target']\n",
    "ax = ptgt.plot(title='COVID-19 R0', x='date',y='target',style='r--', figsize=(10,8))\n",
    "ax.set_yscale(\"symlog\", linthreshy=1)\n",
    "\n",
    "# Rt= Rt0 * pow(2.0,  -(t-t0)/T)\n",
    "Rt0 = 30\n",
    "Th = 7.5\n",
    "Rt1 = Rt0 * pow(2.0, -(55. - 20. ) / Th) \n",
    "rest = pd.DataFrame([[dfs.iloc[20,0],Rt0],[dfs.iloc[55,0],Rt1]])\n",
    "rest.columns = ['date','Rt_estimation']\n",
    "rest.plot(x='date',y='Rt_estimation', ax=ax, style='b:')\n",
    "#\n",
    "for df, label in dflist:\n",
    "    showResult2(ax, df, label)\n",
    "#\n",
    "ax.grid(True)\n",
    "ax.set_ylim(0,)\n",
    "plt.show()\n",
    "fig = ax.get_figure()\n",
    "fig.savefig(\"R0_{}.png\".format(title))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t:2020-03-01 00:00:00, Rt:6.234284221070908\n",
      "t:2020-03-21 00:00:00, Rt:0.9818382401443373\n",
      "t:2020-04-15 00:00:00, Rt:0.09741069095423288\n",
      "t:2020-04-30 00:00:00, Rt:0.02435267273855822\n",
      "t:2020-05-15 00:00:00, Rt:0.006088168184639555\n"
     ]
    }
   ],
   "source": [
    "dfs = dflist[0][0]\n",
    "adv = 17\n",
    "print(\"t:{}, Rt:{}\".format(dfs.iloc[20+adv,0], Rt0 * pow(2.0, -(20.+adv -20. ) / Th) ))\n",
    "adv = 37\n",
    "print(\"t:{}, Rt:{}\".format(dfs.iloc[20+adv,0], Rt0 * pow(2.0, -(20.+adv -20. ) / Th) ))\n",
    "adv = 62\n",
    "dt = pd.Timedelta('1 days')\n",
    "print(\"t:{}, Rt:{}\".format(dfs.iloc[20,0] + adv*dt, Rt0 * pow(2.0, -(20.+adv -20. ) / Th) ))\n",
    "adv = 77\n",
    "dt = pd.Timedelta('1 days')\n",
    "print(\"t:{}, Rt:{}\".format(dfs.iloc[20,0] + adv*dt, Rt0 * pow(2.0, -(20.+adv -20. ) / Th) ))\n",
    "adv = 92\n",
    "dt = pd.Timedelta('1 days')\n",
    "print(\"t:{}, Rt:{}\".format(dfs.iloc[20,0] + adv*dt, Rt0 * pow(2.0, -(20.+adv -20. ) / Th) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}