Figure3-Outbreak analysis by R.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: 'dplyr'\n",
      "\n",
      "The following objects are masked from 'package:stats':\n",
      "\n",
      "    filter, lag\n",
      "\n",
      "The following objects are masked from 'package:base':\n",
      "\n",
      "    intersect, setdiff, setequal, union\n",
      "\n",
      "\n",
      "Attaching package: 'MASS'\n",
      "\n",
      "The following object is masked from 'package:dplyr':\n",
      "\n",
      "    select\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(dplyr)\n",
    "library(MASS)    # for negative binomial regression\n",
    "library(readr)   # read csv file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <div class='info-circle alert alert-block alert-info'> Load combined data</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Parsed with column specification:\n",
      "cols(\n",
      "  .default = col_double(),\n",
      "  city_cn = col_character(),\n",
      "  city = col_character(),\n",
      "  province_cn = col_character(),\n",
      "  province = col_character(),\n",
      "  tier = col_character()\n",
      ")\n",
      "See spec(...) for full column specifications.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>366</li>\n",
       "\t<li>92</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 366\n",
       "\\item 92\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 366\n",
       "2. 92\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 366  92"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reg_df <- read_csv('Data/combined_data.csv')\n",
    "dim(reg_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 6 × 92</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>city_cn</th><th scope=col>city</th><th scope=col>city_id</th><th scope=col>province_cn</th><th scope=col>province</th><th scope=col>tier</th><th scope=col>lat</th><th scope=col>lng</th><th scope=col>population</th><th scope=col>city_tier1</th><th scope=col>...</th><th scope=col>date_0324</th><th scope=col>date_0325</th><th scope=col>date_0326</th><th scope=col>date_0327</th><th scope=col>date_0328</th><th scope=col>date_0329</th><th scope=col>date_0330</th><th scope=col>date_0331</th><th scope=col>date_0401</th><th scope=col>hubei_city</th></tr>\n",
       "\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>...</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>北京市  </td><td>Beijing     </td><td>110000</td><td>北京市</td><td>Beijing</td><td>Tier 1    </td><td>39.91092</td><td>116.4134</td><td>2153.60</td><td>1</td><td>...</td><td>559</td><td>565</td><td>569</td><td>572</td><td>576</td><td>577</td><td>580</td><td>580</td><td>582</td><td>0</td></tr>\n",
       "\t<tr><td>天津市  </td><td>Tianjin     </td><td>120000</td><td>天津市</td><td>Tianjin</td><td>new Tier 1</td><td>39.09367</td><td>117.2095</td><td>1561.83</td><td>1</td><td>...</td><td>145</td><td>147</td><td>152</td><td>156</td><td>163</td><td>166</td><td>174</td><td>174</td><td>176</td><td>0</td></tr>\n",
       "\t<tr><td>石家庄市</td><td>Shijiazhuang</td><td>130100</td><td>河北省</td><td>Hebei  </td><td>Tier 2    </td><td>38.04831</td><td>114.5215</td><td>1039.42</td><td>0</td><td>...</td><td> 29</td><td> 29</td><td> 29</td><td> 29</td><td> 29</td><td> 29</td><td> 29</td><td> 29</td><td> 29</td><td>0</td></tr>\n",
       "\t<tr><td>唐山市  </td><td>Tangshan    </td><td>130200</td><td>河北省</td><td>Hebei  </td><td>Tier 3    </td><td>39.63658</td><td>118.1865</td><td> 796.40</td><td>0</td><td>...</td><td> 58</td><td> 58</td><td> 58</td><td> 58</td><td> 58</td><td> 58</td><td> 58</td><td> 58</td><td> 58</td><td>0</td></tr>\n",
       "\t<tr><td>秦皇岛市</td><td>Qinhuangdao </td><td>130300</td><td>河北省</td><td>Hebei  </td><td>Tier 3    </td><td>39.94175</td><td>119.6085</td><td> 314.63</td><td>0</td><td>...</td><td> 10</td><td> 10</td><td> 10</td><td> 10</td><td> 10</td><td> 10</td><td> 10</td><td> 10</td><td> 10</td><td>0</td></tr>\n",
       "\t<tr><td>邯郸市  </td><td>Handan      </td><td>130400</td><td>河北省</td><td>Hebei  </td><td>Tier 3    </td><td>36.63126</td><td>114.5456</td><td> 954.97</td><td>0</td><td>...</td><td> 32</td><td> 32</td><td> 32</td><td> 32</td><td> 32</td><td> 32</td><td> 32</td><td> 32</td><td> 32</td><td>0</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 6 × 92\n",
       "\\begin{tabular}{r|llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n",
       " city\\_cn & city & city\\_id & province\\_cn & province & tier & lat & lng & population & city\\_tier1 & city\\_tier2 & city\\_tier3 & city\\_tier4 & city\\_tier5 & distance & local\\_flow & pagerank & mean\\_intensity & intensity\\_2019 & intensity\\_2020 & intensity\\_reduction & date\\_0123 & date\\_0124 & date\\_0125 & date\\_0126 & date\\_0127 & date\\_0128 & date\\_0129 & date\\_0130 & date\\_0131 & date\\_0201 & date\\_0202 & date\\_0203 & date\\_0204 & date\\_0205 & date\\_0206 & date\\_0207 & date\\_0208 & date\\_0209 & date\\_0210 & date\\_0211 & date\\_0212 & date\\_0213 & date\\_0214 & date\\_0215 & date\\_0216 & date\\_0217 & date\\_0218 & date\\_0219 & date\\_0220 & date\\_0221 & date\\_0222 & date\\_0223 & date\\_0224 & date\\_0225 & date\\_0226 & date\\_0227 & date\\_0228 & date\\_0229 & date\\_0301 & date\\_0302 & date\\_0303 & date\\_0304 & date\\_0305 & date\\_0306 & date\\_0307 & date\\_0308 & date\\_0309 & date\\_0310 & date\\_0311 & date\\_0312 & date\\_0313 & date\\_0314 & date\\_0315 & date\\_0316 & date\\_0317 & date\\_0318 & date\\_0319 & date\\_0320 & date\\_0321 & date\\_0322 & date\\_0323 & date\\_0324 & date\\_0325 & date\\_0326 & date\\_0327 & date\\_0328 & date\\_0329 & date\\_0330 & date\\_0331 & date\\_0401 & hubei\\_city\\\\\n",
       " <chr> & <chr> & <dbl> & <chr> & <chr> & <chr> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t 北京市   & Beijing      & 110000 & 北京市 & Beijing & Tier 1     & 39.91092 & 116.4134 & 2153.60 & 1 & 0 & 0 & 0 & 0 & 1050.6402 & 184.353501 & 0.015574171 & 5.677559 & 331.9212 & 241.5429 & 90.3783 & 26 & 36 & 51 & 67 & 80 & 96 & 112 & 132 & 156 & 183 & 212 & 228 & 253 & 274 & 297 & 315 & 326 & 337 & 342 & 352 & 366 & 372 & 375 & 380 & 381 & 387 & 393 & 395 & 396 & 399 & 399 & 399 & 400 & 400 & 410 & 410 & 411 & 413 & 414 & 414 & 417 & 418 & 422 & 426 & 428 & 428 & 429 & 435 & 435 & 436 & 437 & 442 & 446 & 455 & 458 & 479 & 485 & 499 & 512 & 522 & 554 & 559 & 565 & 569 & 572 & 576 & 577 & 580 & 580 & 582 & 0\\\\\n",
       "\t 天津市   & Tianjin      & 120000 & 天津市 & Tianjin & new Tier 1 & 39.09367 & 117.2095 & 1561.83 & 1 & 0 & 0 & 0 & 0 &  978.8402 &  30.484535 & 0.005574312 & 5.610923 & 356.6680 & 258.5126 & 98.1554 &  0 &  8 & 10 & 14 & 23 & 24 &  27 &  31 &  32 &  41 &  48 &  60 &  67 &  69 &  79 &  81 &  88 &  91 &  95 & 107 & 113 & 119 & 120 & 122 & 125 & 127 & 128 & 130 & 131 & 133 & 135 & 135 & 135 & 135 & 135 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 136 & 137 & 137 & 137 & 141 & 142 & 145 & 147 & 152 & 156 & 163 & 166 & 174 & 174 & 176 & 0\\\\\n",
       "\t 石家庄市 & Shijiazhuang & 130100 & 河北省 & Hebei   & Tier 2     & 38.04831 & 114.5215 & 1039.42 & 0 & 1 & 0 & 0 & 0 &  826.6277 &  30.892180 & 0.005190617 & 5.479823 & 346.9245 & 260.6993 & 86.2252 &  1 &  4 &  5 &  7 &  9 & 10 &  11 &  11 &  12 &  13 &  16 &  21 &  21 &  24 &  24 &  24 &  24 &  25 &  27 &  27 &  27 &  27 &  27 &  28 &  28 &  28 &  28 &  28 &  28 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 &  29 & 0\\\\\n",
       "\t 唐山市   & Tangshan     & 130200 & 河北省 & Hebei   & Tier 3     & 39.63658 & 118.1865 &  796.40 & 0 & 0 & 1 & 0 & 0 & 1062.8057 &   9.256908 & 0.002692226 & 5.583632 & 355.8822 & 288.3937 & 67.4885 &  0 &  0 &  0 &  0 &  1 &  3 &   3 &   7 &   8 &  11 &  11 &  11 &  12 &  17 &  19 &  23 &  27 &  28 &  31 &  33 &  35 &  43 &  46 &  51 &  51 &  51 &  55 &  56 &  57 &  57 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 &  58 & 0\\\\\n",
       "\t 秦皇岛市 & Qinhuangdao  & 130300 & 河北省 & Hebei   & Tier 3     & 39.94175 & 119.6085 &  314.63 & 0 & 0 & 1 & 0 & 0 & 1142.5763 &   4.230623 & 0.001596582 & 5.384295 & 333.1633 & 280.1461 & 53.0172 &  0 &  0 &  0 &  0 &  0 &  0 &   0 &   0 &   1 &   1 &   2 &   2 &   3 &   3 &   5 &   7 &   7 &   7 &   9 &   9 &   9 &   9 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 &  10 & 0\\\\\n",
       "\t 邯郸市   & Handan       & 130400 & 河北省 & Hebei   & Tier 3     & 36.63126 & 114.5456 &  954.97 & 0 & 0 & 1 & 0 & 0 &  669.4844 &  25.975888 & 0.003952657 & 5.281932 & 354.6887 & 277.3320 & 77.3567 &  0 &  0 &  2 &  2 &  4 &  4 &   6 &   7 &  10 &  11 &  11 &  14 &  14 &  18 &  22 &  26 &  27 &  27 &  27 &  29 &  31 &  31 &  31 &  31 &  31 &  31 &  31 &  31 &  31 &  31 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 &  32 & 0\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 6 × 92\n",
       "\n",
       "| city_cn &lt;chr&gt; | city &lt;chr&gt; | city_id &lt;dbl&gt; | province_cn &lt;chr&gt; | province &lt;chr&gt; | tier &lt;chr&gt; | lat &lt;dbl&gt; | lng &lt;dbl&gt; | population &lt;dbl&gt; | city_tier1 &lt;dbl&gt; | ... ... | date_0324 &lt;dbl&gt; | date_0325 &lt;dbl&gt; | date_0326 &lt;dbl&gt; | date_0327 &lt;dbl&gt; | date_0328 &lt;dbl&gt; | date_0329 &lt;dbl&gt; | date_0330 &lt;dbl&gt; | date_0331 &lt;dbl&gt; | date_0401 &lt;dbl&gt; | hubei_city &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| 北京市   | Beijing      | 110000 | 北京市 | Beijing | Tier 1     | 39.91092 | 116.4134 | 2153.60 | 1 | ... | 559 | 565 | 569 | 572 | 576 | 577 | 580 | 580 | 582 | 0 |\n",
       "| 天津市   | Tianjin      | 120000 | 天津市 | Tianjin | new Tier 1 | 39.09367 | 117.2095 | 1561.83 | 1 | ... | 145 | 147 | 152 | 156 | 163 | 166 | 174 | 174 | 176 | 0 |\n",
       "| 石家庄市 | Shijiazhuang | 130100 | 河北省 | Hebei   | Tier 2     | 38.04831 | 114.5215 | 1039.42 | 0 | ... |  29 |  29 |  29 |  29 |  29 |  29 |  29 |  29 |  29 | 0 |\n",
       "| 唐山市   | Tangshan     | 130200 | 河北省 | Hebei   | Tier 3     | 39.63658 | 118.1865 |  796.40 | 0 | ... |  58 |  58 |  58 |  58 |  58 |  58 |  58 |  58 |  58 | 0 |\n",
       "| 秦皇岛市 | Qinhuangdao  | 130300 | 河北省 | Hebei   | Tier 3     | 39.94175 | 119.6085 |  314.63 | 0 | ... |  10 |  10 |  10 |  10 |  10 |  10 |  10 |  10 |  10 | 0 |\n",
       "| 邯郸市   | Handan       | 130400 | 河北省 | Hebei   | Tier 3     | 36.63126 | 114.5456 |  954.97 | 0 | ... |  32 |  32 |  32 |  32 |  32 |  32 |  32 |  32 |  32 | 0 |\n",
       "\n"
      ],
      "text/plain": [
       "  city_cn  city         city_id province_cn province tier       lat     \n",
       "1 北京市   Beijing      110000  北京市      Beijing  Tier 1     39.91092\n",
       "2 天津市   Tianjin      120000  天津市      Tianjin  new Tier 1 39.09367\n",
       "3 石家庄市 Shijiazhuang 130100  河北省      Hebei    Tier 2     38.04831\n",
       "4 唐山市   Tangshan     130200  河北省      Hebei    Tier 3     39.63658\n",
       "5 秦皇岛市 Qinhuangdao  130300  河北省      Hebei    Tier 3     39.94175\n",
       "6 邯郸市   Handan       130400  河北省      Hebei    Tier 3     36.63126\n",
       "  lng      population city_tier1 ... date_0324 date_0325 date_0326 date_0327\n",
       "1 116.4134 2153.60    1          ... 559       565       569       572      \n",
       "2 117.2095 1561.83    1          ... 145       147       152       156      \n",
       "3 114.5215 1039.42    0          ...  29        29        29        29      \n",
       "4 118.1865  796.40    0          ...  58        58        58        58      \n",
       "5 119.6085  314.63    0          ...  10        10        10        10      \n",
       "6 114.5456  954.97    0          ...  32        32        32        32      \n",
       "  date_0328 date_0329 date_0330 date_0331 date_0401 hubei_city\n",
       "1 576       577       580       580       582       0         \n",
       "2 163       166       174       174       176       0         \n",
       "3  29        29        29        29        29       0         \n",
       "4  58        58        58        58        58       0         \n",
       "5  10        10        10        10        10       0         \n",
       "6  32        32        32        32        32       0         "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head(reg_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter Wuhan\n",
    "reg_df <-\n",
    "reg_df %>% \n",
    "    filter(city != 'Wuhan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <div class=\"alert alert-block alert-danger\">Outbreak Estimates</div>\n",
    "\n",
    "* COVID-19 outbreak in different areas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Outbreak case Prediction using Negative Binomial Model\n",
    "\n",
    "1. For the outbreak in Wuhan, a Binomial regression model is fitted first;\n",
    "2. When the outbreak area changes, the corresponding predictors (e.g., human flow from and the distance to the outbreak city) also vary with the change;\n",
    "3. Assuming that the spatial spread pattern of COVID-19 remains unchanged (which is a very strong assumption);\n",
    "4. Based on the adjusted predictors, estimations are then made based on the fitted model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formula:\n",
    "\n",
    "cumulative_cases ~ exp(log_local_flow + log_population + log_distance + mean_intensity)\n",
    "\n",
    "1. local_flow: human flow from the outbreak city to each city (vary with the outbreak area)\n",
    "2. population: the population of each city (unchanged)\n",
    "3. distance: distance from each city to the outbreak city (vary with the outbreak area)\n",
    "4. mean_intensity: mean intra-city activity intensity (unchanged)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "#### step 1: model fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm.nb(formula = cumulative_cases ~ log_local_flow + log_population + \n",
       "    log_distance + mean_intensity, data = reg_df, init.theta = 2.051967158, \n",
       "    link = log)\n",
       "\n",
       "Deviance Residuals: \n",
       "    Min       1Q   Median       3Q      Max  \n",
       "-2.7898  -0.9566  -0.3890   0.2487   4.4077  \n",
       "\n",
       "Coefficients:\n",
       "               Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)     0.02088    0.83420   0.025   0.9800    \n",
       "log_local_flow  0.79124    0.05517  14.341  < 2e-16 ***\n",
       "log_population  0.32241    0.05741   5.616 1.95e-08 ***\n",
       "log_distance   -0.16561    0.09794  -1.691   0.0908 .  \n",
       "mean_intensity  0.07306    0.06863   1.065   0.2871    \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for Negative Binomial(2.052) family taken to be 1)\n",
       "\n",
       "    Null deviance: 2332.22  on 363  degrees of freedom\n",
       "Residual deviance:  430.82  on 359  degrees of freedom\n",
       "  (1 observation deleted due to missingness)\n",
       "AIC: 2750.7\n",
       "\n",
       "Number of Fisher Scoring iterations: 1\n",
       "\n",
       "\n",
       "              Theta:  2.052 \n",
       "          Std. Err.:  0.192 \n",
       "\n",
       " 2 x log-likelihood:  -2738.712 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# negative Binomial model\n",
    "# in the MASS package\n",
    "\n",
    "# cases on '2020-02-09'\n",
    "reg_df$cumulative_cases <- reg_df$date_0209\n",
    "\n",
    "reg_df <- reg_df %>% \n",
    "                mutate(log_local_flow = log(local_flow + 1)) %>%\n",
    "                mutate(log_population = log(population)) %>%\n",
    "                mutate(log_distance = log(distance + 1))\n",
    "\n",
    "\n",
    "negb_model <- glm.nb(cumulative_cases ~ \n",
    "                    log_local_flow + \n",
    "                    log_population + \n",
    "                    log_distance + \n",
    "                    mean_intensity, \n",
    "              data = reg_df\n",
    "                   )\n",
    "\n",
    "summary(negb_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.487337234451204"
      ],
      "text/latex": [
       "0.487337234451204"
      ],
      "text/markdown": [
       "0.487337234451204"
      ],
      "text/plain": [
       "[1] 0.4873372"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "1/negb_model$theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "#### step 2: model prediciton when the outbreak city changes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### case for Chengdu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A df[,6]: 6 × 6</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>city</th><th scope=col>city_id</th><th scope=col>local_flow</th><th scope=col>population</th><th scope=col>distance</th><th scope=col>mean_intensity</th></tr>\n",
       "\t<tr><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>Beijing     </td><td>110000</td><td>237.216012</td><td>2153.60</td><td>1518.037</td><td>5.677559</td></tr>\n",
       "\t<tr><td>Tianjin     </td><td>120000</td><td> 40.574199</td><td>1561.83</td><td>1519.337</td><td>5.610923</td></tr>\n",
       "\t<tr><td>Shijiazhuang</td><td>130100</td><td> 36.633426</td><td>1039.42</td><td>1261.403</td><td>5.479823</td></tr>\n",
       "\t<tr><td>Tangshan    </td><td>130200</td><td> 11.380362</td><td> 796.40</td><td>1622.828</td><td>5.583632</td></tr>\n",
       "\t<tr><td>Qinhuangdao </td><td>130300</td><td>  5.695457</td><td> 314.63</td><td>1743.796</td><td>5.384295</td></tr>\n",
       "\t<tr><td>Handan      </td><td>130400</td><td> 20.725976</td><td> 954.97</td><td>1174.229</td><td>5.281932</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A df[,6]: 6 × 6\n",
       "\\begin{tabular}{r|llllll}\n",
       " city & city\\_id & local\\_flow & population & distance & mean\\_intensity\\\\\n",
       " <fct> & <int> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t Beijing      & 110000 & 237.216012 & 2153.60 & 1518.037 & 5.677559\\\\\n",
       "\t Tianjin      & 120000 &  40.574199 & 1561.83 & 1519.337 & 5.610923\\\\\n",
       "\t Shijiazhuang & 130100 &  36.633426 & 1039.42 & 1261.403 & 5.479823\\\\\n",
       "\t Tangshan     & 130200 &  11.380362 &  796.40 & 1622.828 & 5.583632\\\\\n",
       "\t Qinhuangdao  & 130300 &   5.695457 &  314.63 & 1743.796 & 5.384295\\\\\n",
       "\t Handan       & 130400 &  20.725976 &  954.97 & 1174.229 & 5.281932\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A df[,6]: 6 × 6\n",
       "\n",
       "| city &lt;fct&gt; | city_id &lt;int&gt; | local_flow &lt;dbl&gt; | population &lt;dbl&gt; | distance &lt;dbl&gt; | mean_intensity &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|\n",
       "| Beijing      | 110000 | 237.216012 | 2153.60 | 1518.037 | 5.677559 |\n",
       "| Tianjin      | 120000 |  40.574199 | 1561.83 | 1519.337 | 5.610923 |\n",
       "| Shijiazhuang | 130100 |  36.633426 | 1039.42 | 1261.403 | 5.479823 |\n",
       "| Tangshan     | 130200 |  11.380362 |  796.40 | 1622.828 | 5.583632 |\n",
       "| Qinhuangdao  | 130300 |   5.695457 |  314.63 | 1743.796 | 5.384295 |\n",
       "| Handan       | 130400 |  20.725976 |  954.97 | 1174.229 | 5.281932 |\n",
       "\n"
      ],
      "text/plain": [
       "  city         city_id local_flow population distance mean_intensity\n",
       "1 Beijing      110000  237.216012 2153.60    1518.037 5.677559      \n",
       "2 Tianjin      120000   40.574199 1561.83    1519.337 5.610923      \n",
       "3 Shijiazhuang 130100   36.633426 1039.42    1261.403 5.479823      \n",
       "4 Tangshan     130200   11.380362  796.40    1622.828 5.583632      \n",
       "5 Qinhuangdao  130300    5.695457  314.63    1743.796 5.384295      \n",
       "6 Handan       130400   20.725976  954.97    1174.229 5.281932      "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# case example\n",
    "city <- 'Chengdu'\n",
    "\n",
    "nb_predict_base_df = read.csv(paste0('Data/outbreak_base_data/nb_base_', city, '_for_prediction.csv'))\n",
    "\n",
    "head(nb_predict_base_df)\n",
    "\n",
    "# \n",
    "# the human flow from Chengdu to Beijing is: 237.216012\n",
    "# the distance from Beijing to Chengdu is: 1518.037"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_predict_base_df <- nb_predict_base_df %>% \n",
    "                mutate(log_local_flow = log(local_flow + 1)) %>%\n",
    "                mutate(log_population = log(population)) %>%\n",
    "                mutate(log_distance = log(distance + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chengdu ,  30492.83 ,  30495 \n"
     ]
    }
   ],
   "source": [
    "predicted_cases <- predict(negb_model, nb_predict_base_df, type = \"response\")\n",
    "\n",
    "# excluding NA: the outbreak city itsself and cities with missing values\n",
    "raw_predicted = sum(predicted_cases, na.rm=TRUE)\n",
    "round_predicted = sum(round(predicted_cases), na.rm=TRUE)\n",
    "\n",
    "cat(city, ', ', raw_predicted, ', ', round_predicted, '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### multiple cities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beijing ,  43722.09 ,  43719 \n",
      "Shanghai ,  40726.77 ,  40737 \n",
      "Tianjin ,  20707.53 ,  20707 \n",
      "Chongqing ,  18023.89 ,  18016 \n",
      "Zhengzhou ,  26258.65 ,  26253 \n",
      "Changsha ,  21162.38 ,  21168 \n",
      "Hefei ,  16151.74 ,  16154 \n",
      "Guangzhou ,  43819.59 ,  43824 \n",
      "Shenzhen ,  40363.96 ,  40363 \n",
      "Nanchang ,  11134.13 ,  11125 \n",
      "Nanjing ,  21501.9 ,  21495 \n",
      "Hangzhou ,  27439.6 ,  27432 \n",
      "Shenyang ,  12126.04 ,  12131 \n",
      "Suzhou ,  33250.31 ,  33251 \n",
      "Xian ,  22668.72 ,  22664 \n",
      "Chengdu ,  30492.83 ,  30495 \n"
     ]
    }
   ],
   "source": [
    "\n",
    "cities <- c('Beijing', 'Shanghai', 'Tianjin', 'Chongqing', \n",
    "            'Zhengzhou', 'Changsha', 'Hefei', 'Guangzhou', \n",
    "            'Shenzhen', 'Nanchang', 'Nanjing', 'Hangzhou', \n",
    "            'Shenyang', 'Suzhou', 'Xian', 'Chengdu')\n",
    "\n",
    "cities_cn <- c('北京市', '上海市', '天津市', '重庆市', \n",
    "              '郑州市', '长沙市', '合肥市', '广州市', \n",
    "              '深圳市', '南昌市', '南京市', '杭州市', \n",
    "              '沈阳市', '苏州市', '西安市', '成都市')\n",
    "\n",
    "raw_predicts <- c()\n",
    "round_predicts <- c()\n",
    "\n",
    "for (city in cities)\n",
    "{\n",
    "    nb_predict_base_df = read.csv(paste0('Data/outbreak_base_data/nb_base_', \n",
    "                                   city, '_for_prediction.csv'))\n",
    "    nb_predict_base_df <- nb_predict_base_df %>% \n",
    "                mutate(log_local_flow = log(local_flow + 1)) %>%\n",
    "                mutate(log_population = log(population)) %>%\n",
    "                mutate(log_distance = log(distance + 1))\n",
    "\n",
    "    predicted_cases <- predict(negb_model, nb_predict_base_df, type = \"response\")\n",
    "    raw_predicted = sum(predicted_cases, na.rm=TRUE)\n",
    "    round_predicted = sum(round(predicted_cases), na.rm=TRUE)\n",
    "    cat(city, ', ', raw_predicted, ', ', round_predicted, '\\n')\n",
    "    # cities <- rbind(cities, city)\n",
    "    raw_predicts <- rbind(raw_predicts, raw_predicted)\n",
    "    round_predicts <- rbind(round_predicts, round_predicted)\n",
    "}\n",
    "\n",
    "# cities_cn\n",
    "# cities\n",
    "# raw_predicts\n",
    "# round_predicts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A df[,4]: 16 × 4</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>city_cn</th><th scope=col>city</th><th scope=col>raw_predict</th><th scope=col>predicted_cumulative_cases</th></tr>\n",
       "\t<tr><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>北京市</td><td>Beijing  </td><td>43722.09</td><td>43719</td></tr>\n",
       "\t<tr><td>上海市</td><td>Shanghai </td><td>40726.77</td><td>40737</td></tr>\n",
       "\t<tr><td>天津市</td><td>Tianjin  </td><td>20707.53</td><td>20707</td></tr>\n",
       "\t<tr><td>重庆市</td><td>Chongqing</td><td>18023.89</td><td>18016</td></tr>\n",
       "\t<tr><td>郑州市</td><td>Zhengzhou</td><td>26258.65</td><td>26253</td></tr>\n",
       "\t<tr><td>长沙市</td><td>Changsha </td><td>21162.38</td><td>21168</td></tr>\n",
       "\t<tr><td>合肥市</td><td>Hefei    </td><td>16151.74</td><td>16154</td></tr>\n",
       "\t<tr><td>广州市</td><td>Guangzhou</td><td>43819.59</td><td>43824</td></tr>\n",
       "\t<tr><td>深圳市</td><td>Shenzhen </td><td>40363.96</td><td>40363</td></tr>\n",
       "\t<tr><td>南昌市</td><td>Nanchang </td><td>11134.13</td><td>11125</td></tr>\n",
       "\t<tr><td>南京市</td><td>Nanjing  </td><td>21501.90</td><td>21495</td></tr>\n",
       "\t<tr><td>杭州市</td><td>Hangzhou </td><td>27439.60</td><td>27432</td></tr>\n",
       "\t<tr><td>沈阳市</td><td>Shenyang </td><td>12126.04</td><td>12131</td></tr>\n",
       "\t<tr><td>苏州市</td><td>Suzhou   </td><td>33250.31</td><td>33251</td></tr>\n",
       "\t<tr><td>西安市</td><td>Xian     </td><td>22668.72</td><td>22664</td></tr>\n",
       "\t<tr><td>成都市</td><td>Chengdu  </td><td>30492.83</td><td>30495</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A df[,4]: 16 × 4\n",
       "\\begin{tabular}{r|llll}\n",
       " city\\_cn & city & raw\\_predict & predicted\\_cumulative\\_cases\\\\\n",
       " <fct> & <fct> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t 北京市 & Beijing   & 43722.09 & 43719\\\\\n",
       "\t 上海市 & Shanghai  & 40726.77 & 40737\\\\\n",
       "\t 天津市 & Tianjin   & 20707.53 & 20707\\\\\n",
       "\t 重庆市 & Chongqing & 18023.89 & 18016\\\\\n",
       "\t 郑州市 & Zhengzhou & 26258.65 & 26253\\\\\n",
       "\t 长沙市 & Changsha  & 21162.38 & 21168\\\\\n",
       "\t 合肥市 & Hefei     & 16151.74 & 16154\\\\\n",
       "\t 广州市 & Guangzhou & 43819.59 & 43824\\\\\n",
       "\t 深圳市 & Shenzhen  & 40363.96 & 40363\\\\\n",
       "\t 南昌市 & Nanchang  & 11134.13 & 11125\\\\\n",
       "\t 南京市 & Nanjing   & 21501.90 & 21495\\\\\n",
       "\t 杭州市 & Hangzhou  & 27439.60 & 27432\\\\\n",
       "\t 沈阳市 & Shenyang  & 12126.04 & 12131\\\\\n",
       "\t 苏州市 & Suzhou    & 33250.31 & 33251\\\\\n",
       "\t 西安市 & Xian      & 22668.72 & 22664\\\\\n",
       "\t 成都市 & Chengdu   & 30492.83 & 30495\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A df[,4]: 16 × 4\n",
       "\n",
       "| city_cn &lt;fct&gt; | city &lt;fct&gt; | raw_predict &lt;dbl&gt; | predicted_cumulative_cases &lt;dbl&gt; |\n",
       "|---|---|---|---|\n",
       "| 北京市 | Beijing   | 43722.09 | 43719 |\n",
       "| 上海市 | Shanghai  | 40726.77 | 40737 |\n",
       "| 天津市 | Tianjin   | 20707.53 | 20707 |\n",
       "| 重庆市 | Chongqing | 18023.89 | 18016 |\n",
       "| 郑州市 | Zhengzhou | 26258.65 | 26253 |\n",
       "| 长沙市 | Changsha  | 21162.38 | 21168 |\n",
       "| 合肥市 | Hefei     | 16151.74 | 16154 |\n",
       "| 广州市 | Guangzhou | 43819.59 | 43824 |\n",
       "| 深圳市 | Shenzhen  | 40363.96 | 40363 |\n",
       "| 南昌市 | Nanchang  | 11134.13 | 11125 |\n",
       "| 南京市 | Nanjing   | 21501.90 | 21495 |\n",
       "| 杭州市 | Hangzhou  | 27439.60 | 27432 |\n",
       "| 沈阳市 | Shenyang  | 12126.04 | 12131 |\n",
       "| 苏州市 | Suzhou    | 33250.31 | 33251 |\n",
       "| 西安市 | Xian      | 22668.72 | 22664 |\n",
       "| 成都市 | Chengdu   | 30492.83 | 30495 |\n",
       "\n"
      ],
      "text/plain": [
       "   city_cn city      raw_predict predicted_cumulative_cases\n",
       "1  北京市  Beijing   43722.09    43719                     \n",
       "2  上海市  Shanghai  40726.77    40737                     \n",
       "3  天津市  Tianjin   20707.53    20707                     \n",
       "4  重庆市  Chongqing 18023.89    18016                     \n",
       "5  郑州市  Zhengzhou 26258.65    26253                     \n",
       "6  长沙市  Changsha  21162.38    21168                     \n",
       "7  合肥市  Hefei     16151.74    16154                     \n",
       "8  广州市  Guangzhou 43819.59    43824                     \n",
       "9  深圳市  Shenzhen  40363.96    40363                     \n",
       "10 南昌市  Nanchang  11134.13    11125                     \n",
       "11 南京市  Nanjing   21501.90    21495                     \n",
       "12 杭州市  Hangzhou  27439.60    27432                     \n",
       "13 沈阳市  Shenyang  12126.04    12131                     \n",
       "14 苏州市  Suzhou    33250.31    33251                     \n",
       "15 西安市  Xian      22668.72    22664                     \n",
       "16 成都市  Chengdu   30492.83    30495                     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "predicted <-\n",
    "data.frame(city_cn=cities_cn, city=cities, \n",
    "           raw_predict=raw_predicts, predicted_cumulative_cases=round_predicts)\n",
    "\n",
    "rownames(predicted) = 1:nrow(predicted)\n",
    "predicted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results to file\n",
    "file_path = file.path(paste0(\"Data/\", \"negb_predicted_outbreak_size_by_area.csv\"))\n",
    "# write.csv(predicted, file = file_path, row.names = FALSE, fileEncoding = \"UTF-8\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}