{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "This is surveillance 1.18.0. For overview type 'help(surveillance)'.\n",
      "\n",
      "Registered S3 method overwritten by 'cli':\n",
      "  method     from    \n",
      "  print.boxx spatstat\n",
      "\n",
      "rstan (Version 2.19.3, GitRev: 2e1f913d3ca3)\n",
      "\n",
      "For execution on a local, multicore CPU with excess RAM we recommend calling\n",
      "options(mc.cores = parallel::detectCores()).\n",
      "To avoid recompilation of unchanged Stan programs, we recommend calling\n",
      "rstan_options(auto_write = TRUE)\n",
      "\n",
      "For improved execution time, we recommend calling\n",
      "Sys.setenv(LOCAL_CPPFLAGS = '-march=corei7 -mtune=corei7')\n",
      "although this causes Stan to throw an error on a few processors.\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "'R version 4.0.0 Patched (2020-05-24 r78561)'"
      ],
      "text/latex": [
       "'R version 4.0.0 Patched (2020-05-24 r78561)'"
      ],
      "text/markdown": [
       "'R version 4.0.0 Patched (2020-05-24 r78561)'"
      ],
      "text/plain": [
       "[1] \"R version 4.0.0 Patched (2020-05-24 r78561)\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "libraries = c(\"dplyr\", \"magrittr\", \"tidyr\", \"surveillance\", \"rstan\") \n",
    "for(x in libraries) { library(x,character.only=TRUE,warn.conflicts=FALSE,quietly=TRUE) }\n",
    "\n",
    "'%&%' = function(x,y)paste0(x,y)\n",
    "\n",
    "R.Version()$version.string\n",
    "\n",
    "options(mc.cores = parallel::detectCores())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "137"
      ],
      "text/latex": [
       "137"
      ],
      "text/markdown": [
       "137"
      ],
      "text/plain": [
       "[1] 137"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A data.frame: 15404 × 6</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>exp_type</th><th scope=col>is_asymptomatic</th><th scope=col>onset</th><th scope=col>time_onset</th><th scope=col>report</th><th scope=col>time_report</th></tr>\n",
       "\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;date&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;date&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-03</td><td> 9</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-14</td><td>20</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-21</td><td>27</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-23</td><td>29</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-22</td><td>28</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-26</td><td>32</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-14</td><td>20</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-20</td><td>26</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-23</td><td>29</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-25</td><td>31</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-24</td><td>30</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-20</td><td>26</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-30</td><td>36</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-01-24</td><td>30</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-24</td><td>30</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-02-08</td><td>45</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-02-02</td><td>39</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-22</td><td>28</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-29</td><td>35</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-31</td><td>37</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-02-03</td><td>40</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-31</td><td>37</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-02-10</td><td>47</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-02-05</td><td>42</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-02-04</td><td>41</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-01-20</td><td>26</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-02-01</td><td>38</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>imported</td><td>0</td><td>2020-02-13</td><td>50</td><td>NA</td><td>NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>1</td><td>NA</td><td>NA</td><td>2020-02-15</td><td>52</td></tr>\n",
       "\t<tr><td>domestic</td><td>1</td><td>NA</td><td>NA</td><td>2020-02-15</td><td>52</td></tr>\n",
       "\t<tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-06</td><td>133</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>1</td><td>NA</td><td> NA</td><td>2020-05-09</td><td>136</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-29</td><td>126</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-06</td><td>133</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-26</td><td>123</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-06</td><td>133</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-03</td><td>130</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-04</td><td>131</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-30</td><td>127</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-07</td><td>134</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-02</td><td>129</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-26</td><td>123</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-24</td><td>121</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-20</td><td>117</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-06</td><td>133</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>1</td><td>NA</td><td> NA</td><td>2020-05-09</td><td>136</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-18</td><td>115</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-02</td><td>129</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-28</td><td>125</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-02</td><td>129</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-06</td><td>133</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-26</td><td>123</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-01</td><td>128</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-06</td><td>133</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-03</td><td>130</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-05</td><td>132</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-04-30</td><td>127</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>0</td><td>2020-05-06</td><td>133</td><td>NA</td><td> NA</td></tr>\n",
       "\t<tr><td>domestic</td><td>1</td><td>NA</td><td> NA</td><td>2020-05-09</td><td>136</td></tr>\n",
       "\t<tr><td>domestic</td><td>1</td><td>NA</td><td> NA</td><td>2020-05-09</td><td>136</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 15404 × 6\n",
       "\\begin{tabular}{llllll}\n",
       " exp\\_type & is\\_asymptomatic & onset & time\\_onset & report & time\\_report\\\\\n",
       " <chr> & <int> & <date> & <dbl> & <date> & <dbl>\\\\\n",
       "\\hline\n",
       "\t imported & 0 & 2020-01-03 &  9 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-14 & 20 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-21 & 27 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-23 & 29 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-22 & 28 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-26 & 32 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-14 & 20 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-20 & 26 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-23 & 29 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-25 & 31 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-24 & 30 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-20 & 26 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-30 & 36 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-01-24 & 30 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-24 & 30 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-02-08 & 45 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-02-02 & 39 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-22 & 28 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-29 & 35 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-31 & 37 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-02-03 & 40 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-31 & 37 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-02-10 & 47 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-02-05 & 42 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-02-04 & 41 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-01-20 & 26 & NA & NA\\\\\n",
       "\t domestic & 0 & 2020-02-01 & 38 & NA & NA\\\\\n",
       "\t imported & 0 & 2020-02-13 & 50 & NA & NA\\\\\n",
       "\t domestic & 1 & NA & NA & 2020-02-15 & 52\\\\\n",
       "\t domestic & 1 & NA & NA & 2020-02-15 & 52\\\\\n",
       "\t ... & ... & ... & ... & ... & ...\\\\\n",
       "\t domestic & 0 & 2020-05-06 & 133 & NA &  NA\\\\\n",
       "\t domestic & 1 & NA &  NA & 2020-05-09 & 136\\\\\n",
       "\t domestic & 0 & 2020-04-29 & 126 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-06 & 133 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-26 & 123 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-06 & 133 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-03 & 130 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-04 & 131 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-30 & 127 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-07 & 134 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-02 & 129 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-26 & 123 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-24 & 121 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-20 & 117 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-06 & 133 & NA &  NA\\\\\n",
       "\t domestic & 1 & NA &  NA & 2020-05-09 & 136\\\\\n",
       "\t domestic & 0 & 2020-04-18 & 115 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-02 & 129 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-28 & 125 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-02 & 129 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-06 & 133 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-26 & 123 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-01 & 128 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-06 & 133 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-03 & 130 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-05 & 132 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-04-30 & 127 & NA &  NA\\\\\n",
       "\t domestic & 0 & 2020-05-06 & 133 & NA &  NA\\\\\n",
       "\t domestic & 1 & NA &  NA & 2020-05-09 & 136\\\\\n",
       "\t domestic & 1 & NA &  NA & 2020-05-09 & 136\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 15404 × 6\n",
       "\n",
       "| exp_type &lt;chr&gt; | is_asymptomatic &lt;int&gt; | onset &lt;date&gt; | time_onset &lt;dbl&gt; | report &lt;date&gt; | time_report &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|\n",
       "| imported | 0 | 2020-01-03 |  9 | NA | NA |\n",
       "| imported | 0 | 2020-01-14 | 20 | NA | NA |\n",
       "| imported | 0 | 2020-01-21 | 27 | NA | NA |\n",
       "| imported | 0 | 2020-01-23 | 29 | NA | NA |\n",
       "| imported | 0 | 2020-01-22 | 28 | NA | NA |\n",
       "| imported | 0 | 2020-01-26 | 32 | NA | NA |\n",
       "| domestic | 0 | 2020-01-14 | 20 | NA | NA |\n",
       "| domestic | 0 | 2020-01-20 | 26 | NA | NA |\n",
       "| domestic | 0 | 2020-01-23 | 29 | NA | NA |\n",
       "| domestic | 0 | 2020-01-25 | 31 | NA | NA |\n",
       "| imported | 0 | 2020-01-24 | 30 | NA | NA |\n",
       "| domestic | 0 | 2020-01-20 | 26 | NA | NA |\n",
       "| imported | 0 | 2020-01-30 | 36 | NA | NA |\n",
       "| imported | 0 | 2020-01-24 | 30 | NA | NA |\n",
       "| domestic | 0 | 2020-01-24 | 30 | NA | NA |\n",
       "| domestic | 0 | 2020-02-08 | 45 | NA | NA |\n",
       "| domestic | 0 | 2020-02-02 | 39 | NA | NA |\n",
       "| domestic | 0 | 2020-01-22 | 28 | NA | NA |\n",
       "| domestic | 0 | 2020-01-29 | 35 | NA | NA |\n",
       "| domestic | 0 | 2020-01-31 | 37 | NA | NA |\n",
       "| imported | 0 | 2020-02-03 | 40 | NA | NA |\n",
       "| domestic | 0 | 2020-01-31 | 37 | NA | NA |\n",
       "| domestic | 0 | 2020-02-10 | 47 | NA | NA |\n",
       "| domestic | 0 | 2020-02-05 | 42 | NA | NA |\n",
       "| domestic | 0 | 2020-02-04 | 41 | NA | NA |\n",
       "| domestic | 0 | 2020-01-20 | 26 | NA | NA |\n",
       "| domestic | 0 | 2020-02-01 | 38 | NA | NA |\n",
       "| imported | 0 | 2020-02-13 | 50 | NA | NA |\n",
       "| domestic | 1 | NA | NA | 2020-02-15 | 52 |\n",
       "| domestic | 1 | NA | NA | 2020-02-15 | 52 |\n",
       "| ... | ... | ... | ... | ... | ... |\n",
       "| domestic | 0 | 2020-05-06 | 133 | NA |  NA |\n",
       "| domestic | 1 | NA |  NA | 2020-05-09 | 136 |\n",
       "| domestic | 0 | 2020-04-29 | 126 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-06 | 133 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-26 | 123 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-06 | 133 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-03 | 130 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-04 | 131 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-30 | 127 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-07 | 134 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-02 | 129 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-26 | 123 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-24 | 121 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-20 | 117 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-06 | 133 | NA |  NA |\n",
       "| domestic | 1 | NA |  NA | 2020-05-09 | 136 |\n",
       "| domestic | 0 | 2020-04-18 | 115 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-02 | 129 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-28 | 125 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-02 | 129 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-06 | 133 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-26 | 123 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-01 | 128 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-06 | 133 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-03 | 130 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-05 | 132 | NA |  NA |\n",
       "| domestic | 0 | 2020-04-30 | 127 | NA |  NA |\n",
       "| domestic | 0 | 2020-05-06 | 133 | NA |  NA |\n",
       "| domestic | 1 | NA |  NA | 2020-05-09 | 136 |\n",
       "| domestic | 1 | NA |  NA | 2020-05-09 | 136 |\n",
       "\n"
      ],
      "text/plain": [
       "      exp_type is_asymptomatic onset      time_onset report     time_report\n",
       "1     imported 0               2020-01-03  9         <NA>       NA         \n",
       "2     imported 0               2020-01-14 20         <NA>       NA         \n",
       "3     imported 0               2020-01-21 27         <NA>       NA         \n",
       "4     imported 0               2020-01-23 29         <NA>       NA         \n",
       "5     imported 0               2020-01-22 28         <NA>       NA         \n",
       "6     imported 0               2020-01-26 32         <NA>       NA         \n",
       "7     domestic 0               2020-01-14 20         <NA>       NA         \n",
       "8     domestic 0               2020-01-20 26         <NA>       NA         \n",
       "9     domestic 0               2020-01-23 29         <NA>       NA         \n",
       "10    domestic 0               2020-01-25 31         <NA>       NA         \n",
       "11    imported 0               2020-01-24 30         <NA>       NA         \n",
       "12    domestic 0               2020-01-20 26         <NA>       NA         \n",
       "13    imported 0               2020-01-30 36         <NA>       NA         \n",
       "14    imported 0               2020-01-24 30         <NA>       NA         \n",
       "15    domestic 0               2020-01-24 30         <NA>       NA         \n",
       "16    domestic 0               2020-02-08 45         <NA>       NA         \n",
       "17    domestic 0               2020-02-02 39         <NA>       NA         \n",
       "18    domestic 0               2020-01-22 28         <NA>       NA         \n",
       "19    domestic 0               2020-01-29 35         <NA>       NA         \n",
       "20    domestic 0               2020-01-31 37         <NA>       NA         \n",
       "21    imported 0               2020-02-03 40         <NA>       NA         \n",
       "22    domestic 0               2020-01-31 37         <NA>       NA         \n",
       "23    domestic 0               2020-02-10 47         <NA>       NA         \n",
       "24    domestic 0               2020-02-05 42         <NA>       NA         \n",
       "25    domestic 0               2020-02-04 41         <NA>       NA         \n",
       "26    domestic 0               2020-01-20 26         <NA>       NA         \n",
       "27    domestic 0               2020-02-01 38         <NA>       NA         \n",
       "28    imported 0               2020-02-13 50         <NA>       NA         \n",
       "29    domestic 1               <NA>       NA         2020-02-15 52         \n",
       "30    domestic 1               <NA>       NA         2020-02-15 52         \n",
       "...   ...      ...             ...        ...        ...        ...        \n",
       "15375 domestic 0               2020-05-06 133        <NA>        NA        \n",
       "15376 domestic 1               <NA>        NA        2020-05-09 136        \n",
       "15377 domestic 0               2020-04-29 126        <NA>        NA        \n",
       "15378 domestic 0               2020-05-06 133        <NA>        NA        \n",
       "15379 domestic 0               2020-04-26 123        <NA>        NA        \n",
       "15380 domestic 0               2020-05-06 133        <NA>        NA        \n",
       "15381 domestic 0               2020-05-03 130        <NA>        NA        \n",
       "15382 domestic 0               2020-05-04 131        <NA>        NA        \n",
       "15383 domestic 0               2020-04-30 127        <NA>        NA        \n",
       "15384 domestic 0               2020-05-07 134        <NA>        NA        \n",
       "15385 domestic 0               2020-05-02 129        <NA>        NA        \n",
       "15386 domestic 0               2020-04-26 123        <NA>        NA        \n",
       "15387 domestic 0               2020-04-24 121        <NA>        NA        \n",
       "15388 domestic 0               2020-04-20 117        <NA>        NA        \n",
       "15389 domestic 0               2020-05-06 133        <NA>        NA        \n",
       "15390 domestic 1               <NA>        NA        2020-05-09 136        \n",
       "15391 domestic 0               2020-04-18 115        <NA>        NA        \n",
       "15392 domestic 0               2020-05-02 129        <NA>        NA        \n",
       "15393 domestic 0               2020-04-28 125        <NA>        NA        \n",
       "15394 domestic 0               2020-05-02 129        <NA>        NA        \n",
       "15395 domestic 0               2020-05-06 133        <NA>        NA        \n",
       "15396 domestic 0               2020-04-26 123        <NA>        NA        \n",
       "15397 domestic 0               2020-05-01 128        <NA>        NA        \n",
       "15398 domestic 0               2020-05-06 133        <NA>        NA        \n",
       "15399 domestic 0               2020-05-03 130        <NA>        NA        \n",
       "15400 domestic 0               2020-05-05 132        <NA>        NA        \n",
       "15401 domestic 0               2020-04-30 127        <NA>        NA        \n",
       "15402 domestic 0               2020-05-06 133        <NA>        NA        \n",
       "15403 domestic 1               <NA>        NA        2020-05-09 136        \n",
       "15404 domestic 1               <NA>        NA        2020-05-09 136        "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "datestar = as.Date(\"2020-05-10\")\n",
    "datemin = as.Date(\"2019-12-25\") # particular choice\n",
    "\n",
    "read.csv(\"../data/JapaneseDataCOVID19 (\"%&%format(datestar,\"%y%m%d\")%&%\").csv\") %>%\n",
    "    mutate(report = if_else(is.na(confirmed), reported, confirmed), \n",
    "           # report and onset as Dates\n",
    "           report = as.Date(report), \n",
    "           onset = as.Date(onset),\n",
    "           # time_report and time_onset as the number of days since datemin\n",
    "           time_onset = as.numeric(onset - datemin),  \n",
    "           time_report = as.numeric(report - datemin)) %>%\n",
    "    # removing original column as no need anymore\n",
    "    select(-confirmed, -reported) -> df\n",
    "\n",
    "# nicer order ot the columns\n",
    "df %<>% select(exp_type, is_asymptomatic, onset, time_onset, report, time_report)\n",
    "\n",
    "# number of days between the cutoff time and datemin\n",
    "(tstar = as.numeric(datestar - datemin))\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`summarise()` regrouping output by 'exp_type' (override with `.groups` argument)\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 6 × 3</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>time</th><th scope=col>domestic</th><th scope=col>imported</th></tr>\n",
       "\t<tr><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td> 22</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td> 23</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td> 24</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td> 25</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td> 34</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>136</td><td>0</td><td>0</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 6 × 3\n",
       "\\begin{tabular}{lll}\n",
       " time & domestic & imported\\\\\n",
       " <dbl> & <int> & <int>\\\\\n",
       "\\hline\n",
       "\t  22 & 0 & 0\\\\\n",
       "\t  23 & 0 & 0\\\\\n",
       "\t  24 & 0 & 0\\\\\n",
       "\t  25 & 0 & 0\\\\\n",
       "\t  34 & 0 & 0\\\\\n",
       "\t 136 & 0 & 0\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 6 × 3\n",
       "\n",
       "| time &lt;dbl&gt; | domestic &lt;int&gt; | imported &lt;int&gt; |\n",
       "|---|---|---|\n",
       "|  22 | 0 | 0 |\n",
       "|  23 | 0 | 0 |\n",
       "|  24 | 0 | 0 |\n",
       "|  25 | 0 | 0 |\n",
       "|  34 | 0 | 0 |\n",
       "| 136 | 0 | 0 |\n",
       "\n"
      ],
      "text/plain": [
       "  time domestic imported\n",
       "1  22  0        0       \n",
       "2  23  0        0       \n",
       "3  24  0        0       \n",
       "4  25  0        0       \n",
       "5  34  0        0       \n",
       "6 136  0        0       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# At the moment we have a line-list, but we also want to have the number of imported and domestic cases per day\n",
    "# so we need to count and change the dataframe from long to wide format \n",
    "df %>%\n",
    "    # respectively to the day of illness onset (time_onset)\n",
    "    select(exp_type, time_onset) %>% \n",
    "    # to count the cases we group by exp_type (domestic vs imported) and by time_onset\n",
    "    group_by(exp_type, time_onset) %>% \n",
    "    summarize(n = n()) %>%\n",
    "    # then we transform the dataframe to the wide format\n",
    "    spread(exp_type, n) %>% \n",
    "    # however some days are missing because there were no cases on those days (e.g. at the beginning of the outbreak)\n",
    "    # so we account for all subsequent days from datemin to datestar\n",
    "    right_join(expand_grid(time_onset = seq(from = 1, to = tstar-1, by = 1)), by='time_onset') %>% \n",
    "    # those missing days would have NAs, so we replace them with zeros\n",
    "    replace(is.na(.), 0) %>%\n",
    "    # for simplicity renaming the time columns\n",
    "    rename(time = time_onset) -> Df_cases\n",
    "\n",
    "Df_cases %>% tail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`summarise()` regrouping output by 'exp_type' (override with `.groups` argument)\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 6 × 5</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>time</th><th scope=col>domestic_unobs</th><th scope=col>imported_unobs</th><th scope=col>domestic</th><th scope=col>imported</th></tr>\n",
       "\t<tr><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>131</td><td>27</td><td>0</td><td>28</td><td>0</td></tr>\n",
       "\t<tr><td>132</td><td>27</td><td>0</td><td>19</td><td>0</td></tr>\n",
       "\t<tr><td>133</td><td>15</td><td>0</td><td>17</td><td>0</td></tr>\n",
       "\t<tr><td>134</td><td>11</td><td>0</td><td> 7</td><td>0</td></tr>\n",
       "\t<tr><td>135</td><td>12</td><td>0</td><td> 4</td><td>0</td></tr>\n",
       "\t<tr><td>136</td><td>11</td><td>0</td><td> 0</td><td>0</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 6 × 5\n",
       "\\begin{tabular}{lllll}\n",
       " time & domestic\\_unobs & imported\\_unobs & domestic & imported\\\\\n",
       " <dbl> & <int> & <int> & <int> & <int>\\\\\n",
       "\\hline\n",
       "\t 131 & 27 & 0 & 28 & 0\\\\\n",
       "\t 132 & 27 & 0 & 19 & 0\\\\\n",
       "\t 133 & 15 & 0 & 17 & 0\\\\\n",
       "\t 134 & 11 & 0 &  7 & 0\\\\\n",
       "\t 135 & 12 & 0 &  4 & 0\\\\\n",
       "\t 136 & 11 & 0 &  0 & 0\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 6 × 5\n",
       "\n",
       "| time &lt;dbl&gt; | domestic_unobs &lt;int&gt; | imported_unobs &lt;int&gt; | domestic &lt;int&gt; | imported &lt;int&gt; |\n",
       "|---|---|---|---|---|\n",
       "| 131 | 27 | 0 | 28 | 0 |\n",
       "| 132 | 27 | 0 | 19 | 0 |\n",
       "| 133 | 15 | 0 | 17 | 0 |\n",
       "| 134 | 11 | 0 |  7 | 0 |\n",
       "| 135 | 12 | 0 |  4 | 0 |\n",
       "| 136 | 11 | 0 |  0 | 0 |\n",
       "\n"
      ],
      "text/plain": [
       "  time domestic_unobs imported_unobs domestic imported\n",
       "1 131  27             0              28       0       \n",
       "2 132  27             0              19       0       \n",
       "3 133  15             0              17       0       \n",
       "4 134  11             0               7       0       \n",
       "5 135  12             0               4       0       \n",
       "6 136  11             0               0       0       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# because some cases do not have date of illness onset, we still want to count them and probably do the backprojection later\n",
    "# so we repeat the previous procedure but by the day of report\n",
    "df %>% \n",
    "    filter(is.na(time_onset)) %>% \n",
    "    select(exp_type, time_report) %>% \n",
    "    group_by(exp_type, time_report) %>% \n",
    "    summarize(n = n()) %>% \n",
    "    spread(exp_type, n) %>% \n",
    "    # we mark those cases with \"_unobs\"-postfix\n",
    "    rename(imported_unobs = imported, \n",
    "           domestic_unobs = domestic, \n",
    "           time = time_report) %>% \n",
    "    # joining with the previously obtained dataframe\n",
    "    right_join(Df_cases, by='time') %>% \n",
    "    # again some days have zero counts, so they will be NAs after the right_join => making them zeros\n",
    "    replace(is.na(.), 0) %>%\n",
    "    arrange(time) -> Df_cases\n",
    "\n",
    "Df_cases %>% tail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"purple\">Backprojection</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 6 × 5</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>time</th><th scope=col>domestic_unobs</th><th scope=col>imported_unobs</th><th scope=col>domestic</th><th scope=col>imported</th></tr>\n",
       "\t<tr><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>142</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>143</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>144</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>145</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>146</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>147</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 6 × 5\n",
       "\\begin{tabular}{lllll}\n",
       " time & domestic\\_unobs & imported\\_unobs & domestic & imported\\\\\n",
       " <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t 142 & 0 & 0 & 0 & 0\\\\\n",
       "\t 143 & 0 & 0 & 0 & 0\\\\\n",
       "\t 144 & 0 & 0 & 0 & 0\\\\\n",
       "\t 145 & 0 & 0 & 0 & 0\\\\\n",
       "\t 146 & 0 & 0 & 0 & 0\\\\\n",
       "\t 147 & 0 & 0 & 0 & 0\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 6 × 5\n",
       "\n",
       "| time &lt;dbl&gt; | domestic_unobs &lt;dbl&gt; | imported_unobs &lt;dbl&gt; | domestic &lt;dbl&gt; | imported &lt;dbl&gt; |\n",
       "|---|---|---|---|---|\n",
       "| 142 | 0 | 0 | 0 | 0 |\n",
       "| 143 | 0 | 0 | 0 | 0 |\n",
       "| 144 | 0 | 0 | 0 | 0 |\n",
       "| 145 | 0 | 0 | 0 | 0 |\n",
       "| 146 | 0 | 0 | 0 | 0 |\n",
       "| 147 | 0 | 0 | 0 | 0 |\n",
       "\n"
      ],
      "text/plain": [
       "  time domestic_unobs imported_unobs domestic imported\n",
       "1 142  0              0              0        0       \n",
       "2 143  0              0              0        0       \n",
       "3 144  0              0              0        0       \n",
       "4 145  0              0              0        0       \n",
       "5 146  0              0              0        0       \n",
       "6 147  0              0              0        0       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# if the observed counts are close to the right end, the backprojection algorithm can be unstable for them \n",
    "# as outlined by the developer of the package\n",
    "# so we will play a bit safer and put additional 10 days after datestar filled with zeros\n",
    "df_cases = Df_cases %>% rbind(data.frame(time=seq(tstar,tstar+10,1), imported=0, imported_unobs=0, domestic=0, domestic_unobs=0))\n",
    "df_cases %>% tail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 30 × 7</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>time</th><th scope=col>domestic_unobs</th><th scope=col>imported_unobs</th><th scope=col>domestic</th><th scope=col>imported</th><th scope=col>imported_backproj</th><th scope=col>domestic_backproj</th></tr>\n",
       "\t<tr><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>118</td><td>115</td><td>0</td><td>155</td><td>0</td><td>0</td><td>98.108752345</td></tr>\n",
       "\t<tr><td>119</td><td> 71</td><td>0</td><td>162</td><td>0</td><td>0</td><td>93.958232329</td></tr>\n",
       "\t<tr><td>120</td><td> 83</td><td>0</td><td>106</td><td>0</td><td>0</td><td>86.820185261</td></tr>\n",
       "\t<tr><td>121</td><td> 66</td><td>0</td><td>147</td><td>0</td><td>0</td><td>76.029308817</td></tr>\n",
       "\t<tr><td>122</td><td> 47</td><td>0</td><td>115</td><td>0</td><td>0</td><td>62.129701932</td></tr>\n",
       "\t<tr><td>123</td><td> 22</td><td>0</td><td>104</td><td>0</td><td>0</td><td>47.032092685</td></tr>\n",
       "\t<tr><td>124</td><td> 47</td><td>0</td><td>121</td><td>0</td><td>0</td><td>32.748434689</td></tr>\n",
       "\t<tr><td>125</td><td>117</td><td>0</td><td> 94</td><td>0</td><td>0</td><td>20.166454186</td></tr>\n",
       "\t<tr><td>126</td><td> 49</td><td>0</td><td> 69</td><td>0</td><td>0</td><td>10.333683929</td></tr>\n",
       "\t<tr><td>127</td><td> 44</td><td>0</td><td> 77</td><td>0</td><td>0</td><td> 4.207846904</td></tr>\n",
       "\t<tr><td>128</td><td> 42</td><td>0</td><td> 70</td><td>0</td><td>0</td><td> 1.331383489</td></tr>\n",
       "\t<tr><td>129</td><td> 44</td><td>0</td><td> 59</td><td>0</td><td>0</td><td> 0.322522824</td></tr>\n",
       "\t<tr><td>130</td><td> 37</td><td>0</td><td> 39</td><td>0</td><td>0</td><td> 0.057996779</td></tr>\n",
       "\t<tr><td>131</td><td> 27</td><td>0</td><td> 28</td><td>0</td><td>0</td><td> 0.007152493</td></tr>\n",
       "\t<tr><td>132</td><td> 27</td><td>0</td><td> 19</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>133</td><td> 15</td><td>0</td><td> 17</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>134</td><td> 11</td><td>0</td><td>  7</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>135</td><td> 12</td><td>0</td><td>  4</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>136</td><td> 11</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>137</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>138</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>139</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>140</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>141</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>142</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>143</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>144</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>145</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>146</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>147</td><td>  0</td><td>0</td><td>  0</td><td>0</td><td>0</td><td> 0.000000000</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 30 × 7\n",
       "\\begin{tabular}{lllllll}\n",
       " time & domestic\\_unobs & imported\\_unobs & domestic & imported & imported\\_backproj & domestic\\_backproj\\\\\n",
       " <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t 118 & 115 & 0 & 155 & 0 & 0 & 98.108752345\\\\\n",
       "\t 119 &  71 & 0 & 162 & 0 & 0 & 93.958232329\\\\\n",
       "\t 120 &  83 & 0 & 106 & 0 & 0 & 86.820185261\\\\\n",
       "\t 121 &  66 & 0 & 147 & 0 & 0 & 76.029308817\\\\\n",
       "\t 122 &  47 & 0 & 115 & 0 & 0 & 62.129701932\\\\\n",
       "\t 123 &  22 & 0 & 104 & 0 & 0 & 47.032092685\\\\\n",
       "\t 124 &  47 & 0 & 121 & 0 & 0 & 32.748434689\\\\\n",
       "\t 125 & 117 & 0 &  94 & 0 & 0 & 20.166454186\\\\\n",
       "\t 126 &  49 & 0 &  69 & 0 & 0 & 10.333683929\\\\\n",
       "\t 127 &  44 & 0 &  77 & 0 & 0 &  4.207846904\\\\\n",
       "\t 128 &  42 & 0 &  70 & 0 & 0 &  1.331383489\\\\\n",
       "\t 129 &  44 & 0 &  59 & 0 & 0 &  0.322522824\\\\\n",
       "\t 130 &  37 & 0 &  39 & 0 & 0 &  0.057996779\\\\\n",
       "\t 131 &  27 & 0 &  28 & 0 & 0 &  0.007152493\\\\\n",
       "\t 132 &  27 & 0 &  19 & 0 & 0 &  0.000000000\\\\\n",
       "\t 133 &  15 & 0 &  17 & 0 & 0 &  0.000000000\\\\\n",
       "\t 134 &  11 & 0 &   7 & 0 & 0 &  0.000000000\\\\\n",
       "\t 135 &  12 & 0 &   4 & 0 & 0 &  0.000000000\\\\\n",
       "\t 136 &  11 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 137 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 138 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 139 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 140 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 141 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 142 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 143 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 144 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 145 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 146 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\t 147 &   0 & 0 &   0 & 0 & 0 &  0.000000000\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 30 × 7\n",
       "\n",
       "| time &lt;dbl&gt; | domestic_unobs &lt;dbl&gt; | imported_unobs &lt;dbl&gt; | domestic &lt;dbl&gt; | imported &lt;dbl&gt; | imported_backproj &lt;dbl&gt; | domestic_backproj &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|\n",
       "| 118 | 115 | 0 | 155 | 0 | 0 | 98.108752345 |\n",
       "| 119 |  71 | 0 | 162 | 0 | 0 | 93.958232329 |\n",
       "| 120 |  83 | 0 | 106 | 0 | 0 | 86.820185261 |\n",
       "| 121 |  66 | 0 | 147 | 0 | 0 | 76.029308817 |\n",
       "| 122 |  47 | 0 | 115 | 0 | 0 | 62.129701932 |\n",
       "| 123 |  22 | 0 | 104 | 0 | 0 | 47.032092685 |\n",
       "| 124 |  47 | 0 | 121 | 0 | 0 | 32.748434689 |\n",
       "| 125 | 117 | 0 |  94 | 0 | 0 | 20.166454186 |\n",
       "| 126 |  49 | 0 |  69 | 0 | 0 | 10.333683929 |\n",
       "| 127 |  44 | 0 |  77 | 0 | 0 |  4.207846904 |\n",
       "| 128 |  42 | 0 |  70 | 0 | 0 |  1.331383489 |\n",
       "| 129 |  44 | 0 |  59 | 0 | 0 |  0.322522824 |\n",
       "| 130 |  37 | 0 |  39 | 0 | 0 |  0.057996779 |\n",
       "| 131 |  27 | 0 |  28 | 0 | 0 |  0.007152493 |\n",
       "| 132 |  27 | 0 |  19 | 0 | 0 |  0.000000000 |\n",
       "| 133 |  15 | 0 |  17 | 0 | 0 |  0.000000000 |\n",
       "| 134 |  11 | 0 |   7 | 0 | 0 |  0.000000000 |\n",
       "| 135 |  12 | 0 |   4 | 0 | 0 |  0.000000000 |\n",
       "| 136 |  11 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 137 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 138 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 139 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 140 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 141 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 142 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 143 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 144 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 145 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 146 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "| 147 |   0 | 0 |   0 | 0 | 0 |  0.000000000 |\n",
       "\n"
      ],
      "text/plain": [
       "   time domestic_unobs imported_unobs domestic imported imported_backproj\n",
       "1  118  115            0              155      0        0                \n",
       "2  119   71            0              162      0        0                \n",
       "3  120   83            0              106      0        0                \n",
       "4  121   66            0              147      0        0                \n",
       "5  122   47            0              115      0        0                \n",
       "6  123   22            0              104      0        0                \n",
       "7  124   47            0              121      0        0                \n",
       "8  125  117            0               94      0        0                \n",
       "9  126   49            0               69      0        0                \n",
       "10 127   44            0               77      0        0                \n",
       "11 128   42            0               70      0        0                \n",
       "12 129   44            0               59      0        0                \n",
       "13 130   37            0               39      0        0                \n",
       "14 131   27            0               28      0        0                \n",
       "15 132   27            0               19      0        0                \n",
       "16 133   15            0               17      0        0                \n",
       "17 134   11            0                7      0        0                \n",
       "18 135   12            0                4      0        0                \n",
       "19 136   11            0                0      0        0                \n",
       "20 137    0            0                0      0        0                \n",
       "21 138    0            0                0      0        0                \n",
       "22 139    0            0                0      0        0                \n",
       "23 140    0            0                0      0        0                \n",
       "24 141    0            0                0      0        0                \n",
       "25 142    0            0                0      0        0                \n",
       "26 143    0            0                0      0        0                \n",
       "27 144    0            0                0      0        0                \n",
       "28 145    0            0                0      0        0                \n",
       "29 146    0            0                0      0        0                \n",
       "30 147    0            0                0      0        0                \n",
       "   domestic_backproj\n",
       "1  98.108752345     \n",
       "2  93.958232329     \n",
       "3  86.820185261     \n",
       "4  76.029308817     \n",
       "5  62.129701932     \n",
       "6  47.032092685     \n",
       "7  32.748434689     \n",
       "8  20.166454186     \n",
       "9  10.333683929     \n",
       "10  4.207846904     \n",
       "11  1.331383489     \n",
       "12  0.322522824     \n",
       "13  0.057996779     \n",
       "14  0.007152493     \n",
       "15  0.000000000     \n",
       "16  0.000000000     \n",
       "17  0.000000000     \n",
       "18  0.000000000     \n",
       "19  0.000000000     \n",
       "20  0.000000000     \n",
       "21  0.000000000     \n",
       "22  0.000000000     \n",
       "23  0.000000000     \n",
       "24  0.000000000     \n",
       "25  0.000000000     \n",
       "26  0.000000000     \n",
       "27  0.000000000     \n",
       "28  0.000000000     \n",
       "29  0.000000000     \n",
       "30  0.000000000     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "K = nrow(df_cases)\n",
    "\n",
    "# Incubation time period is adopted from [Linton et al 2020]\n",
    "# other research groups give quite similar values of the mean 5.2 days\n",
    "# see https://github.com/aakhmetz/COVID19IncubationPeriod\n",
    "inc_fit = list(meanlog=1.519, sdlog=0.615)\n",
    "incubation_probability = plnorm(1:K, inc_fit$meanlog, inc_fit$sdlog) - plnorm(1:K-1, inc_fit$meanlog, inc_fit$sdlog)\n",
    "# the surveillance packages requires to keep zero at the beginning\n",
    "inc_pmf = c(0,incubation_probability[1:41])\n",
    "# smoothing parameter values for backprojection of imported and domestic\n",
    "k_used = c(2,2)\n",
    "\n",
    "# if the backprojected value falls below this threshold, we set it to zero \n",
    "# in this case we will avoid arbitrarily small values which make no sense\n",
    "epsilon = 1e-3\n",
    "\n",
    "# backprojection of imported cases\n",
    "# parameters of the backprojection follows the parameters adopted in the tutorial for the package\n",
    "# the backprojection is done only if there are imported cases in the database\n",
    "if (sum(df_cases$imported)>0) {\n",
    "    sts = new(\"sts\", epoch=df_cases$time, observed=df_cases$imported)\n",
    "    bpnp.control = list(k = k_used[1], eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                        B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, \n",
    "                        eq3a.method = c(\"R\",\"C\"))\n",
    "    sts_bp = backprojNP(sts, incu.pmf=inc_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "    # output result\n",
    "    df_cases['imported_backproj'] = upperbound(sts_bp) %>% as.numeric %>% replace(. < epsilon, 0) \n",
    "    df_cases['imported_backproj'] = df_cases['imported_backproj']/sum(df_cases['imported_backproj'])*sum(df_cases$imported)\n",
    "} else \n",
    "    df_cases['imported_backproj'] = 0\n",
    "\n",
    "## backprojection of domestic cases (=domesically acquired infections)\n",
    "if (sum(df_cases$domestic)>0) {\n",
    "    sts = new(\"sts\", epoch=df_cases$time, observed=df_cases$domestic)\n",
    "    bpnp.control = list(k = k_used[2], eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                        B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, \n",
    "                        eq3a.method = c(\"R\",\"C\"))\n",
    "    sts_bp = backprojNP(sts, incu.pmf=inc_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "    df_cases['domestic_backproj'] = upperbound(sts_bp) %>% as.numeric %>% replace(. < epsilon, 0) \n",
    "    df_cases['domestic_backproj'] = df_cases['domestic_backproj']/sum(df_cases['domestic_backproj'])*sum(df_cases$domestic)\n",
    "} else\n",
    "    df_cases['domestic_backproj'] = 0\n",
    "\n",
    "df_cases %>% tail(30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"purple\">Backprojection including cases with unobserved date of illness onset</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 20 × 9</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>time</th><th scope=col>domestic_unobs</th><th scope=col>imported_unobs</th><th scope=col>domestic</th><th scope=col>imported</th><th scope=col>imported_backproj</th><th scope=col>domestic_backproj</th><th scope=col>imported_backproj_unobs</th><th scope=col>domestic_backproj_unobs</th></tr>\n",
       "\t<tr><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>117</td><td> 98</td><td>0</td><td>219</td><td>0</td><td>0</td><td>1.021782e+02</td><td>0</td><td>55.275342682</td></tr>\n",
       "\t<tr><td>118</td><td>115</td><td>0</td><td>155</td><td>0</td><td>0</td><td>9.810875e+01</td><td>0</td><td>53.085239052</td></tr>\n",
       "\t<tr><td>119</td><td> 71</td><td>0</td><td>162</td><td>0</td><td>0</td><td>9.395823e+01</td><td>0</td><td>51.242101691</td></tr>\n",
       "\t<tr><td>120</td><td> 83</td><td>0</td><td>106</td><td>0</td><td>0</td><td>8.682019e+01</td><td>0</td><td>48.486515636</td></tr>\n",
       "\t<tr><td>121</td><td> 66</td><td>0</td><td>147</td><td>0</td><td>0</td><td>7.602931e+01</td><td>0</td><td>43.398273012</td></tr>\n",
       "\t<tr><td>122</td><td> 47</td><td>0</td><td>115</td><td>0</td><td>0</td><td>6.212970e+01</td><td>0</td><td>35.194230769</td></tr>\n",
       "\t<tr><td>123</td><td> 22</td><td>0</td><td>104</td><td>0</td><td>0</td><td>4.703209e+01</td><td>0</td><td>24.678663446</td></tr>\n",
       "\t<tr><td>124</td><td> 47</td><td>0</td><td>121</td><td>0</td><td>0</td><td>3.274843e+01</td><td>0</td><td>14.449643806</td></tr>\n",
       "\t<tr><td>125</td><td>117</td><td>0</td><td> 94</td><td>0</td><td>0</td><td>2.016645e+01</td><td>0</td><td> 7.017536818</td></tr>\n",
       "\t<tr><td>126</td><td> 49</td><td>0</td><td> 69</td><td>0</td><td>0</td><td>1.033368e+01</td><td>0</td><td> 2.880890504</td></tr>\n",
       "\t<tr><td>127</td><td> 44</td><td>0</td><td> 77</td><td>0</td><td>0</td><td>4.207847e+00</td><td>0</td><td> 1.007775298</td></tr>\n",
       "\t<tr><td>128</td><td> 42</td><td>0</td><td> 70</td><td>0</td><td>0</td><td>1.331383e+00</td><td>0</td><td> 0.293235407</td></tr>\n",
       "\t<tr><td>129</td><td> 44</td><td>0</td><td> 59</td><td>0</td><td>0</td><td>3.225228e-01</td><td>0</td><td> 0.069160341</td></tr>\n",
       "\t<tr><td>130</td><td> 37</td><td>0</td><td> 39</td><td>0</td><td>0</td><td>5.799678e-02</td><td>0</td><td> 0.012909911</td></tr>\n",
       "\t<tr><td>131</td><td> 27</td><td>0</td><td> 28</td><td>0</td><td>0</td><td>7.152493e-03</td><td>0</td><td> 0.001863148</td></tr>\n",
       "\t<tr><td>132</td><td> 27</td><td>0</td><td> 19</td><td>0</td><td>0</td><td>0.000000e+00</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>133</td><td> 15</td><td>0</td><td> 17</td><td>0</td><td>0</td><td>0.000000e+00</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>134</td><td> 11</td><td>0</td><td>  7</td><td>0</td><td>0</td><td>0.000000e+00</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>135</td><td> 12</td><td>0</td><td>  4</td><td>0</td><td>0</td><td>0.000000e+00</td><td>0</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>136</td><td> 11</td><td>0</td><td>  0</td><td>0</td><td>0</td><td>0.000000e+00</td><td>0</td><td> 0.000000000</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 20 × 9\n",
       "\\begin{tabular}{lllllllll}\n",
       " time & domestic\\_unobs & imported\\_unobs & domestic & imported & imported\\_backproj & domestic\\_backproj & imported\\_backproj\\_unobs & domestic\\_backproj\\_unobs\\\\\n",
       " <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t 117 &  98 & 0 & 219 & 0 & 0 & 1.021782e+02 & 0 & 55.275342682\\\\\n",
       "\t 118 & 115 & 0 & 155 & 0 & 0 & 9.810875e+01 & 0 & 53.085239052\\\\\n",
       "\t 119 &  71 & 0 & 162 & 0 & 0 & 9.395823e+01 & 0 & 51.242101691\\\\\n",
       "\t 120 &  83 & 0 & 106 & 0 & 0 & 8.682019e+01 & 0 & 48.486515636\\\\\n",
       "\t 121 &  66 & 0 & 147 & 0 & 0 & 7.602931e+01 & 0 & 43.398273012\\\\\n",
       "\t 122 &  47 & 0 & 115 & 0 & 0 & 6.212970e+01 & 0 & 35.194230769\\\\\n",
       "\t 123 &  22 & 0 & 104 & 0 & 0 & 4.703209e+01 & 0 & 24.678663446\\\\\n",
       "\t 124 &  47 & 0 & 121 & 0 & 0 & 3.274843e+01 & 0 & 14.449643806\\\\\n",
       "\t 125 & 117 & 0 &  94 & 0 & 0 & 2.016645e+01 & 0 &  7.017536818\\\\\n",
       "\t 126 &  49 & 0 &  69 & 0 & 0 & 1.033368e+01 & 0 &  2.880890504\\\\\n",
       "\t 127 &  44 & 0 &  77 & 0 & 0 & 4.207847e+00 & 0 &  1.007775298\\\\\n",
       "\t 128 &  42 & 0 &  70 & 0 & 0 & 1.331383e+00 & 0 &  0.293235407\\\\\n",
       "\t 129 &  44 & 0 &  59 & 0 & 0 & 3.225228e-01 & 0 &  0.069160341\\\\\n",
       "\t 130 &  37 & 0 &  39 & 0 & 0 & 5.799678e-02 & 0 &  0.012909911\\\\\n",
       "\t 131 &  27 & 0 &  28 & 0 & 0 & 7.152493e-03 & 0 &  0.001863148\\\\\n",
       "\t 132 &  27 & 0 &  19 & 0 & 0 & 0.000000e+00 & 0 &  0.000000000\\\\\n",
       "\t 133 &  15 & 0 &  17 & 0 & 0 & 0.000000e+00 & 0 &  0.000000000\\\\\n",
       "\t 134 &  11 & 0 &   7 & 0 & 0 & 0.000000e+00 & 0 &  0.000000000\\\\\n",
       "\t 135 &  12 & 0 &   4 & 0 & 0 & 0.000000e+00 & 0 &  0.000000000\\\\\n",
       "\t 136 &  11 & 0 &   0 & 0 & 0 & 0.000000e+00 & 0 &  0.000000000\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 20 × 9\n",
       "\n",
       "| time &lt;dbl&gt; | domestic_unobs &lt;dbl&gt; | imported_unobs &lt;dbl&gt; | domestic &lt;dbl&gt; | imported &lt;dbl&gt; | imported_backproj &lt;dbl&gt; | domestic_backproj &lt;dbl&gt; | imported_backproj_unobs &lt;dbl&gt; | domestic_backproj_unobs &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|\n",
       "| 117 |  98 | 0 | 219 | 0 | 0 | 1.021782e+02 | 0 | 55.275342682 |\n",
       "| 118 | 115 | 0 | 155 | 0 | 0 | 9.810875e+01 | 0 | 53.085239052 |\n",
       "| 119 |  71 | 0 | 162 | 0 | 0 | 9.395823e+01 | 0 | 51.242101691 |\n",
       "| 120 |  83 | 0 | 106 | 0 | 0 | 8.682019e+01 | 0 | 48.486515636 |\n",
       "| 121 |  66 | 0 | 147 | 0 | 0 | 7.602931e+01 | 0 | 43.398273012 |\n",
       "| 122 |  47 | 0 | 115 | 0 | 0 | 6.212970e+01 | 0 | 35.194230769 |\n",
       "| 123 |  22 | 0 | 104 | 0 | 0 | 4.703209e+01 | 0 | 24.678663446 |\n",
       "| 124 |  47 | 0 | 121 | 0 | 0 | 3.274843e+01 | 0 | 14.449643806 |\n",
       "| 125 | 117 | 0 |  94 | 0 | 0 | 2.016645e+01 | 0 |  7.017536818 |\n",
       "| 126 |  49 | 0 |  69 | 0 | 0 | 1.033368e+01 | 0 |  2.880890504 |\n",
       "| 127 |  44 | 0 |  77 | 0 | 0 | 4.207847e+00 | 0 |  1.007775298 |\n",
       "| 128 |  42 | 0 |  70 | 0 | 0 | 1.331383e+00 | 0 |  0.293235407 |\n",
       "| 129 |  44 | 0 |  59 | 0 | 0 | 3.225228e-01 | 0 |  0.069160341 |\n",
       "| 130 |  37 | 0 |  39 | 0 | 0 | 5.799678e-02 | 0 |  0.012909911 |\n",
       "| 131 |  27 | 0 |  28 | 0 | 0 | 7.152493e-03 | 0 |  0.001863148 |\n",
       "| 132 |  27 | 0 |  19 | 0 | 0 | 0.000000e+00 | 0 |  0.000000000 |\n",
       "| 133 |  15 | 0 |  17 | 0 | 0 | 0.000000e+00 | 0 |  0.000000000 |\n",
       "| 134 |  11 | 0 |   7 | 0 | 0 | 0.000000e+00 | 0 |  0.000000000 |\n",
       "| 135 |  12 | 0 |   4 | 0 | 0 | 0.000000e+00 | 0 |  0.000000000 |\n",
       "| 136 |  11 | 0 |   0 | 0 | 0 | 0.000000e+00 | 0 |  0.000000000 |\n",
       "\n"
      ],
      "text/plain": [
       "   time domestic_unobs imported_unobs domestic imported imported_backproj\n",
       "1  117   98            0              219      0        0                \n",
       "2  118  115            0              155      0        0                \n",
       "3  119   71            0              162      0        0                \n",
       "4  120   83            0              106      0        0                \n",
       "5  121   66            0              147      0        0                \n",
       "6  122   47            0              115      0        0                \n",
       "7  123   22            0              104      0        0                \n",
       "8  124   47            0              121      0        0                \n",
       "9  125  117            0               94      0        0                \n",
       "10 126   49            0               69      0        0                \n",
       "11 127   44            0               77      0        0                \n",
       "12 128   42            0               70      0        0                \n",
       "13 129   44            0               59      0        0                \n",
       "14 130   37            0               39      0        0                \n",
       "15 131   27            0               28      0        0                \n",
       "16 132   27            0               19      0        0                \n",
       "17 133   15            0               17      0        0                \n",
       "18 134   11            0                7      0        0                \n",
       "19 135   12            0                4      0        0                \n",
       "20 136   11            0                0      0        0                \n",
       "   domestic_backproj imported_backproj_unobs domestic_backproj_unobs\n",
       "1  1.021782e+02      0                       55.275342682           \n",
       "2  9.810875e+01      0                       53.085239052           \n",
       "3  9.395823e+01      0                       51.242101691           \n",
       "4  8.682019e+01      0                       48.486515636           \n",
       "5  7.602931e+01      0                       43.398273012           \n",
       "6  6.212970e+01      0                       35.194230769           \n",
       "7  4.703209e+01      0                       24.678663446           \n",
       "8  3.274843e+01      0                       14.449643806           \n",
       "9  2.016645e+01      0                        7.017536818           \n",
       "10 1.033368e+01      0                        2.880890504           \n",
       "11 4.207847e+00      0                        1.007775298           \n",
       "12 1.331383e+00      0                        0.293235407           \n",
       "13 3.225228e-01      0                        0.069160341           \n",
       "14 5.799678e-02      0                        0.012909911           \n",
       "15 7.152493e-03      0                        0.001863148           \n",
       "16 0.000000e+00      0                        0.000000000           \n",
       "17 0.000000e+00      0                        0.000000000           \n",
       "18 0.000000e+00      0                        0.000000000           \n",
       "19 0.000000e+00      0                        0.000000000           \n",
       "20 0.000000e+00      0                        0.000000000           "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Backprojection of cases with unknown (=unobserved) dates of illness onset is done in two subsequent steps\n",
    "# First we backproject to the suggesting date of illness onset\n",
    "# Second we sum up with known observed counts and backproject them together by the incubation period\n",
    "# This is done like that because the surveillance package actually prohibits the heterogeneity in the data\n",
    "\n",
    "# Reporting delay: from onset to confirmation\n",
    "# the values are obtained separetely and here it is assumed to be known\n",
    "onset2report_fit = list(param1 = 1.741, param2 = 8.573)\n",
    "onset2report_probability = pweibull(1:K, onset2report_fit$param1, scale=onset2report_fit$param2) - pweibull(1:K-1, onset2report_fit$param1, scale=onset2report_fit$param2)\n",
    "onset2report_pmf = c(0,onset2report_probability[1:41])\n",
    "\n",
    "## import\n",
    "if (sum(df_cases$imported_unobs)>0) {\n",
    "    sts = new(\"sts\", epoch=df_cases$time, observed=df_cases$imported_unobs)\n",
    "    bpnp.control = list(k = k_used[1], eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                        B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, \n",
    "                        eq3a.method = c(\"R\",\"C\"))\n",
    "    sts_bp = backprojNP(sts, incu.pmf=onset2report_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "    # output result\n",
    "    df_cases['imported_backproj_unobs'] = upperbound(sts_bp) %>% as.numeric %>% replace(. < epsilon, 0) \n",
    "    df_cases['imported_backproj_unobs'] = df_cases['imported_backproj_unobs']/sum(df_cases['imported_backproj_unobs'])*sum(df_cases$imported_unobs)\n",
    "} else \n",
    "    df_cases['imported_backproj_unobs'] = 0\n",
    "\n",
    "## secondary\n",
    "if (sum(df_cases$domestic_unobs)>0) {\n",
    "    sts = new(\"sts\", epoch=df_cases$time, observed=df_cases$domestic_unobs)\n",
    "    bpnp.control = list(k = k_used[2], eps = rep(1e-3,2), iter.max=rep(1000,2), \n",
    "                        Tmark = tstar - 1, #attention\n",
    "                        B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, \n",
    "                        eq3a.method = c(\"C\"))\n",
    "    sts_bp = backprojNP(sts, incu.pmf=onset2report_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "    df_cases['domestic_backproj_unobs'] = upperbound(sts_bp) %>% as.numeric %>% replace(. < epsilon, 0) \n",
    "    df_cases['domestic_backproj_unobs'] = df_cases['domestic_backproj_unobs']/sum(df_cases['domestic_backproj_unobs'])*sum(df_cases$domestic_unobs)\n",
    "} else\n",
    "    df_cases['domestic_backproj_unobs'] = 0\n",
    "df_cases %<>% ungroup\n",
    "\n",
    "df_cases %>% filter(time<tstar) %>% tail(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 30 × 11</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>time</th><th scope=col>domestic_unobs</th><th scope=col>imported_unobs</th><th scope=col>domestic</th><th scope=col>imported</th><th scope=col>imported_backproj</th><th scope=col>domestic_backproj</th><th scope=col>imported_backproj_unobs</th><th scope=col>domestic_backproj_unobs</th><th scope=col>imported_backproj_incl_unobs</th><th scope=col>domestic_backproj_incl_unobs</th></tr>\n",
       "\t<tr><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><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>107</td><td>151</td><td>0</td><td>371</td><td>2</td><td>0.004898637</td><td>2.537108e+02</td><td>0.009760285</td><td>1.009866e+02</td><td>0.003141247</td><td>3.056866e+02</td></tr>\n",
       "\t<tr><td>108</td><td>154</td><td>0</td><td>356</td><td>0</td><td>0.000000000</td><td>2.311702e+02</td><td>0.001660786</td><td>8.799885e+01</td><td>0.000000000</td><td>2.859795e+02</td></tr>\n",
       "\t<tr><td>109</td><td> 56</td><td>0</td><td>283</td><td>2</td><td>0.000000000</td><td>2.061040e+02</td><td>0.000000000</td><td>8.128916e+01</td><td>0.000000000</td><td>2.625939e+02</td></tr>\n",
       "\t<tr><td>110</td><td> 56</td><td>2</td><td>323</td><td>0</td><td>0.000000000</td><td>1.831526e+02</td><td>0.000000000</td><td>7.647833e+01</td><td>0.000000000</td><td>2.398983e+02</td></tr>\n",
       "\t<tr><td>111</td><td>147</td><td>0</td><td>292</td><td>0</td><td>0.000000000</td><td>1.660050e+02</td><td>0.000000000</td><td>7.202435e+01</td><td>0.000000000</td><td>2.220498e+02</td></tr>\n",
       "\t<tr><td>112</td><td>101</td><td>0</td><td>280</td><td>0</td><td>0.000000000</td><td>1.539337e+02</td><td>0.000000000</td><td>6.843435e+01</td><td>0.000000000</td><td>2.087914e+02</td></tr>\n",
       "\t<tr><td>113</td><td> 98</td><td>0</td><td>226</td><td>0</td><td>0.000000000</td><td>1.440922e+02</td><td>0.000000000</td><td>6.556801e+01</td><td>0.000000000</td><td>1.973702e+02</td></tr>\n",
       "\t<tr><td>114</td><td> 91</td><td>0</td><td>225</td><td>0</td><td>0.000000000</td><td>1.338058e+02</td><td>0.000000000</td><td>6.296449e+01</td><td>0.000000000</td><td>1.847542e+02</td></tr>\n",
       "\t<tr><td>115</td><td>119</td><td>0</td><td>216</td><td>0</td><td>0.000000000</td><td>1.216966e+02</td><td>0.000000000</td><td>6.050123e+01</td><td>0.000000000</td><td>1.687925e+02</td></tr>\n",
       "\t<tr><td>116</td><td> 59</td><td>0</td><td>156</td><td>0</td><td>0.000000000</td><td>1.099459e+02</td><td>0.000000000</td><td>5.793726e+01</td><td>0.000000000</td><td>1.508628e+02</td></tr>\n",
       "\t<tr><td>117</td><td> 98</td><td>0</td><td>219</td><td>0</td><td>0.000000000</td><td>1.021782e+02</td><td>0.000000000</td><td>5.527534e+01</td><td>0.000000000</td><td>1.341985e+02</td></tr>\n",
       "\t<tr><td>118</td><td>115</td><td>0</td><td>155</td><td>0</td><td>0.000000000</td><td>9.810875e+01</td><td>0.000000000</td><td>5.308524e+01</td><td>0.000000000</td><td>1.189053e+02</td></tr>\n",
       "\t<tr><td>119</td><td> 71</td><td>0</td><td>162</td><td>0</td><td>0.000000000</td><td>9.395823e+01</td><td>0.000000000</td><td>5.124210e+01</td><td>0.000000000</td><td>1.028381e+02</td></tr>\n",
       "\t<tr><td>120</td><td> 83</td><td>0</td><td>106</td><td>0</td><td>0.000000000</td><td>8.682019e+01</td><td>0.000000000</td><td>4.848652e+01</td><td>0.000000000</td><td>8.550350e+01</td></tr>\n",
       "\t<tr><td>121</td><td> 66</td><td>0</td><td>147</td><td>0</td><td>0.000000000</td><td>7.602931e+01</td><td>0.000000000</td><td>4.339827e+01</td><td>0.000000000</td><td>6.818233e+01</td></tr>\n",
       "\t<tr><td>122</td><td> 47</td><td>0</td><td>115</td><td>0</td><td>0.000000000</td><td>6.212970e+01</td><td>0.000000000</td><td>3.519423e+01</td><td>0.000000000</td><td>5.200660e+01</td></tr>\n",
       "\t<tr><td>123</td><td> 22</td><td>0</td><td>104</td><td>0</td><td>0.000000000</td><td>4.703209e+01</td><td>0.000000000</td><td>2.467866e+01</td><td>0.000000000</td><td>3.783090e+01</td></tr>\n",
       "\t<tr><td>124</td><td> 47</td><td>0</td><td>121</td><td>0</td><td>0.000000000</td><td>3.274843e+01</td><td>0.000000000</td><td>1.444964e+01</td><td>0.000000000</td><td>2.593969e+01</td></tr>\n",
       "\t<tr><td>125</td><td>117</td><td>0</td><td> 94</td><td>0</td><td>0.000000000</td><td>2.016645e+01</td><td>0.000000000</td><td>7.017537e+00</td><td>0.000000000</td><td>1.598117e+01</td></tr>\n",
       "\t<tr><td>126</td><td> 49</td><td>0</td><td> 69</td><td>0</td><td>0.000000000</td><td>1.033368e+01</td><td>0.000000000</td><td>2.880891e+00</td><td>0.000000000</td><td>8.262887e+00</td></tr>\n",
       "\t<tr><td>127</td><td> 44</td><td>0</td><td> 77</td><td>0</td><td>0.000000000</td><td>4.207847e+00</td><td>0.000000000</td><td>1.007775e+00</td><td>0.000000000</td><td>3.407351e+00</td></tr>\n",
       "\t<tr><td>128</td><td> 42</td><td>0</td><td> 70</td><td>0</td><td>0.000000000</td><td>1.331383e+00</td><td>0.000000000</td><td>2.932354e-01</td><td>0.000000000</td><td>1.092231e+00</td></tr>\n",
       "\t<tr><td>129</td><td> 44</td><td>0</td><td> 59</td><td>0</td><td>0.000000000</td><td>3.225228e-01</td><td>0.000000000</td><td>6.916034e-02</td><td>0.000000000</td><td>2.676521e-01</td></tr>\n",
       "\t<tr><td>130</td><td> 37</td><td>0</td><td> 39</td><td>0</td><td>0.000000000</td><td>5.799678e-02</td><td>0.000000000</td><td>1.290991e-02</td><td>0.000000000</td><td>4.857183e-02</td></tr>\n",
       "\t<tr><td>131</td><td> 27</td><td>0</td><td> 28</td><td>0</td><td>0.000000000</td><td>7.152493e-03</td><td>0.000000000</td><td>1.863148e-03</td><td>0.000000000</td><td>6.030619e-03</td></tr>\n",
       "\t<tr><td>132</td><td> 27</td><td>0</td><td> 19</td><td>0</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td></tr>\n",
       "\t<tr><td>133</td><td> 15</td><td>0</td><td> 17</td><td>0</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td></tr>\n",
       "\t<tr><td>134</td><td> 11</td><td>0</td><td>  7</td><td>0</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td></tr>\n",
       "\t<tr><td>135</td><td> 12</td><td>0</td><td>  4</td><td>0</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td></tr>\n",
       "\t<tr><td>136</td><td> 11</td><td>0</td><td>  0</td><td>0</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td><td>0.000000000</td><td>0.000000e+00</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 30 × 11\n",
       "\\begin{tabular}{lllllllllll}\n",
       " time & domestic\\_unobs & imported\\_unobs & domestic & imported & imported\\_backproj & domestic\\_backproj & imported\\_backproj\\_unobs & domestic\\_backproj\\_unobs & imported\\_backproj\\_incl\\_unobs & domestic\\_backproj\\_incl\\_unobs\\\\\n",
       " <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t 107 & 151 & 0 & 371 & 2 & 0.004898637 & 2.537108e+02 & 0.009760285 & 1.009866e+02 & 0.003141247 & 3.056866e+02\\\\\n",
       "\t 108 & 154 & 0 & 356 & 0 & 0.000000000 & 2.311702e+02 & 0.001660786 & 8.799885e+01 & 0.000000000 & 2.859795e+02\\\\\n",
       "\t 109 &  56 & 0 & 283 & 2 & 0.000000000 & 2.061040e+02 & 0.000000000 & 8.128916e+01 & 0.000000000 & 2.625939e+02\\\\\n",
       "\t 110 &  56 & 2 & 323 & 0 & 0.000000000 & 1.831526e+02 & 0.000000000 & 7.647833e+01 & 0.000000000 & 2.398983e+02\\\\\n",
       "\t 111 & 147 & 0 & 292 & 0 & 0.000000000 & 1.660050e+02 & 0.000000000 & 7.202435e+01 & 0.000000000 & 2.220498e+02\\\\\n",
       "\t 112 & 101 & 0 & 280 & 0 & 0.000000000 & 1.539337e+02 & 0.000000000 & 6.843435e+01 & 0.000000000 & 2.087914e+02\\\\\n",
       "\t 113 &  98 & 0 & 226 & 0 & 0.000000000 & 1.440922e+02 & 0.000000000 & 6.556801e+01 & 0.000000000 & 1.973702e+02\\\\\n",
       "\t 114 &  91 & 0 & 225 & 0 & 0.000000000 & 1.338058e+02 & 0.000000000 & 6.296449e+01 & 0.000000000 & 1.847542e+02\\\\\n",
       "\t 115 & 119 & 0 & 216 & 0 & 0.000000000 & 1.216966e+02 & 0.000000000 & 6.050123e+01 & 0.000000000 & 1.687925e+02\\\\\n",
       "\t 116 &  59 & 0 & 156 & 0 & 0.000000000 & 1.099459e+02 & 0.000000000 & 5.793726e+01 & 0.000000000 & 1.508628e+02\\\\\n",
       "\t 117 &  98 & 0 & 219 & 0 & 0.000000000 & 1.021782e+02 & 0.000000000 & 5.527534e+01 & 0.000000000 & 1.341985e+02\\\\\n",
       "\t 118 & 115 & 0 & 155 & 0 & 0.000000000 & 9.810875e+01 & 0.000000000 & 5.308524e+01 & 0.000000000 & 1.189053e+02\\\\\n",
       "\t 119 &  71 & 0 & 162 & 0 & 0.000000000 & 9.395823e+01 & 0.000000000 & 5.124210e+01 & 0.000000000 & 1.028381e+02\\\\\n",
       "\t 120 &  83 & 0 & 106 & 0 & 0.000000000 & 8.682019e+01 & 0.000000000 & 4.848652e+01 & 0.000000000 & 8.550350e+01\\\\\n",
       "\t 121 &  66 & 0 & 147 & 0 & 0.000000000 & 7.602931e+01 & 0.000000000 & 4.339827e+01 & 0.000000000 & 6.818233e+01\\\\\n",
       "\t 122 &  47 & 0 & 115 & 0 & 0.000000000 & 6.212970e+01 & 0.000000000 & 3.519423e+01 & 0.000000000 & 5.200660e+01\\\\\n",
       "\t 123 &  22 & 0 & 104 & 0 & 0.000000000 & 4.703209e+01 & 0.000000000 & 2.467866e+01 & 0.000000000 & 3.783090e+01\\\\\n",
       "\t 124 &  47 & 0 & 121 & 0 & 0.000000000 & 3.274843e+01 & 0.000000000 & 1.444964e+01 & 0.000000000 & 2.593969e+01\\\\\n",
       "\t 125 & 117 & 0 &  94 & 0 & 0.000000000 & 2.016645e+01 & 0.000000000 & 7.017537e+00 & 0.000000000 & 1.598117e+01\\\\\n",
       "\t 126 &  49 & 0 &  69 & 0 & 0.000000000 & 1.033368e+01 & 0.000000000 & 2.880891e+00 & 0.000000000 & 8.262887e+00\\\\\n",
       "\t 127 &  44 & 0 &  77 & 0 & 0.000000000 & 4.207847e+00 & 0.000000000 & 1.007775e+00 & 0.000000000 & 3.407351e+00\\\\\n",
       "\t 128 &  42 & 0 &  70 & 0 & 0.000000000 & 1.331383e+00 & 0.000000000 & 2.932354e-01 & 0.000000000 & 1.092231e+00\\\\\n",
       "\t 129 &  44 & 0 &  59 & 0 & 0.000000000 & 3.225228e-01 & 0.000000000 & 6.916034e-02 & 0.000000000 & 2.676521e-01\\\\\n",
       "\t 130 &  37 & 0 &  39 & 0 & 0.000000000 & 5.799678e-02 & 0.000000000 & 1.290991e-02 & 0.000000000 & 4.857183e-02\\\\\n",
       "\t 131 &  27 & 0 &  28 & 0 & 0.000000000 & 7.152493e-03 & 0.000000000 & 1.863148e-03 & 0.000000000 & 6.030619e-03\\\\\n",
       "\t 132 &  27 & 0 &  19 & 0 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00\\\\\n",
       "\t 133 &  15 & 0 &  17 & 0 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00\\\\\n",
       "\t 134 &  11 & 0 &   7 & 0 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00\\\\\n",
       "\t 135 &  12 & 0 &   4 & 0 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00\\\\\n",
       "\t 136 &  11 & 0 &   0 & 0 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00 & 0.000000000 & 0.000000e+00\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 30 × 11\n",
       "\n",
       "| time &lt;dbl&gt; | domestic_unobs &lt;dbl&gt; | imported_unobs &lt;dbl&gt; | domestic &lt;dbl&gt; | imported &lt;dbl&gt; | imported_backproj &lt;dbl&gt; | domestic_backproj &lt;dbl&gt; | imported_backproj_unobs &lt;dbl&gt; | domestic_backproj_unobs &lt;dbl&gt; | imported_backproj_incl_unobs &lt;dbl&gt; | domestic_backproj_incl_unobs &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| 107 | 151 | 0 | 371 | 2 | 0.004898637 | 2.537108e+02 | 0.009760285 | 1.009866e+02 | 0.003141247 | 3.056866e+02 |\n",
       "| 108 | 154 | 0 | 356 | 0 | 0.000000000 | 2.311702e+02 | 0.001660786 | 8.799885e+01 | 0.000000000 | 2.859795e+02 |\n",
       "| 109 |  56 | 0 | 283 | 2 | 0.000000000 | 2.061040e+02 | 0.000000000 | 8.128916e+01 | 0.000000000 | 2.625939e+02 |\n",
       "| 110 |  56 | 2 | 323 | 0 | 0.000000000 | 1.831526e+02 | 0.000000000 | 7.647833e+01 | 0.000000000 | 2.398983e+02 |\n",
       "| 111 | 147 | 0 | 292 | 0 | 0.000000000 | 1.660050e+02 | 0.000000000 | 7.202435e+01 | 0.000000000 | 2.220498e+02 |\n",
       "| 112 | 101 | 0 | 280 | 0 | 0.000000000 | 1.539337e+02 | 0.000000000 | 6.843435e+01 | 0.000000000 | 2.087914e+02 |\n",
       "| 113 |  98 | 0 | 226 | 0 | 0.000000000 | 1.440922e+02 | 0.000000000 | 6.556801e+01 | 0.000000000 | 1.973702e+02 |\n",
       "| 114 |  91 | 0 | 225 | 0 | 0.000000000 | 1.338058e+02 | 0.000000000 | 6.296449e+01 | 0.000000000 | 1.847542e+02 |\n",
       "| 115 | 119 | 0 | 216 | 0 | 0.000000000 | 1.216966e+02 | 0.000000000 | 6.050123e+01 | 0.000000000 | 1.687925e+02 |\n",
       "| 116 |  59 | 0 | 156 | 0 | 0.000000000 | 1.099459e+02 | 0.000000000 | 5.793726e+01 | 0.000000000 | 1.508628e+02 |\n",
       "| 117 |  98 | 0 | 219 | 0 | 0.000000000 | 1.021782e+02 | 0.000000000 | 5.527534e+01 | 0.000000000 | 1.341985e+02 |\n",
       "| 118 | 115 | 0 | 155 | 0 | 0.000000000 | 9.810875e+01 | 0.000000000 | 5.308524e+01 | 0.000000000 | 1.189053e+02 |\n",
       "| 119 |  71 | 0 | 162 | 0 | 0.000000000 | 9.395823e+01 | 0.000000000 | 5.124210e+01 | 0.000000000 | 1.028381e+02 |\n",
       "| 120 |  83 | 0 | 106 | 0 | 0.000000000 | 8.682019e+01 | 0.000000000 | 4.848652e+01 | 0.000000000 | 8.550350e+01 |\n",
       "| 121 |  66 | 0 | 147 | 0 | 0.000000000 | 7.602931e+01 | 0.000000000 | 4.339827e+01 | 0.000000000 | 6.818233e+01 |\n",
       "| 122 |  47 | 0 | 115 | 0 | 0.000000000 | 6.212970e+01 | 0.000000000 | 3.519423e+01 | 0.000000000 | 5.200660e+01 |\n",
       "| 123 |  22 | 0 | 104 | 0 | 0.000000000 | 4.703209e+01 | 0.000000000 | 2.467866e+01 | 0.000000000 | 3.783090e+01 |\n",
       "| 124 |  47 | 0 | 121 | 0 | 0.000000000 | 3.274843e+01 | 0.000000000 | 1.444964e+01 | 0.000000000 | 2.593969e+01 |\n",
       "| 125 | 117 | 0 |  94 | 0 | 0.000000000 | 2.016645e+01 | 0.000000000 | 7.017537e+00 | 0.000000000 | 1.598117e+01 |\n",
       "| 126 |  49 | 0 |  69 | 0 | 0.000000000 | 1.033368e+01 | 0.000000000 | 2.880891e+00 | 0.000000000 | 8.262887e+00 |\n",
       "| 127 |  44 | 0 |  77 | 0 | 0.000000000 | 4.207847e+00 | 0.000000000 | 1.007775e+00 | 0.000000000 | 3.407351e+00 |\n",
       "| 128 |  42 | 0 |  70 | 0 | 0.000000000 | 1.331383e+00 | 0.000000000 | 2.932354e-01 | 0.000000000 | 1.092231e+00 |\n",
       "| 129 |  44 | 0 |  59 | 0 | 0.000000000 | 3.225228e-01 | 0.000000000 | 6.916034e-02 | 0.000000000 | 2.676521e-01 |\n",
       "| 130 |  37 | 0 |  39 | 0 | 0.000000000 | 5.799678e-02 | 0.000000000 | 1.290991e-02 | 0.000000000 | 4.857183e-02 |\n",
       "| 131 |  27 | 0 |  28 | 0 | 0.000000000 | 7.152493e-03 | 0.000000000 | 1.863148e-03 | 0.000000000 | 6.030619e-03 |\n",
       "| 132 |  27 | 0 |  19 | 0 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 |\n",
       "| 133 |  15 | 0 |  17 | 0 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 |\n",
       "| 134 |  11 | 0 |   7 | 0 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 |\n",
       "| 135 |  12 | 0 |   4 | 0 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 |\n",
       "| 136 |  11 | 0 |   0 | 0 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 | 0.000000000 | 0.000000e+00 |\n",
       "\n"
      ],
      "text/plain": [
       "   time domestic_unobs imported_unobs domestic imported imported_backproj\n",
       "1  107  151            0              371      2        0.004898637      \n",
       "2  108  154            0              356      0        0.000000000      \n",
       "3  109   56            0              283      2        0.000000000      \n",
       "4  110   56            2              323      0        0.000000000      \n",
       "5  111  147            0              292      0        0.000000000      \n",
       "6  112  101            0              280      0        0.000000000      \n",
       "7  113   98            0              226      0        0.000000000      \n",
       "8  114   91            0              225      0        0.000000000      \n",
       "9  115  119            0              216      0        0.000000000      \n",
       "10 116   59            0              156      0        0.000000000      \n",
       "11 117   98            0              219      0        0.000000000      \n",
       "12 118  115            0              155      0        0.000000000      \n",
       "13 119   71            0              162      0        0.000000000      \n",
       "14 120   83            0              106      0        0.000000000      \n",
       "15 121   66            0              147      0        0.000000000      \n",
       "16 122   47            0              115      0        0.000000000      \n",
       "17 123   22            0              104      0        0.000000000      \n",
       "18 124   47            0              121      0        0.000000000      \n",
       "19 125  117            0               94      0        0.000000000      \n",
       "20 126   49            0               69      0        0.000000000      \n",
       "21 127   44            0               77      0        0.000000000      \n",
       "22 128   42            0               70      0        0.000000000      \n",
       "23 129   44            0               59      0        0.000000000      \n",
       "24 130   37            0               39      0        0.000000000      \n",
       "25 131   27            0               28      0        0.000000000      \n",
       "26 132   27            0               19      0        0.000000000      \n",
       "27 133   15            0               17      0        0.000000000      \n",
       "28 134   11            0                7      0        0.000000000      \n",
       "29 135   12            0                4      0        0.000000000      \n",
       "30 136   11            0                0      0        0.000000000      \n",
       "   domestic_backproj imported_backproj_unobs domestic_backproj_unobs\n",
       "1  2.537108e+02      0.009760285             1.009866e+02           \n",
       "2  2.311702e+02      0.001660786             8.799885e+01           \n",
       "3  2.061040e+02      0.000000000             8.128916e+01           \n",
       "4  1.831526e+02      0.000000000             7.647833e+01           \n",
       "5  1.660050e+02      0.000000000             7.202435e+01           \n",
       "6  1.539337e+02      0.000000000             6.843435e+01           \n",
       "7  1.440922e+02      0.000000000             6.556801e+01           \n",
       "8  1.338058e+02      0.000000000             6.296449e+01           \n",
       "9  1.216966e+02      0.000000000             6.050123e+01           \n",
       "10 1.099459e+02      0.000000000             5.793726e+01           \n",
       "11 1.021782e+02      0.000000000             5.527534e+01           \n",
       "12 9.810875e+01      0.000000000             5.308524e+01           \n",
       "13 9.395823e+01      0.000000000             5.124210e+01           \n",
       "14 8.682019e+01      0.000000000             4.848652e+01           \n",
       "15 7.602931e+01      0.000000000             4.339827e+01           \n",
       "16 6.212970e+01      0.000000000             3.519423e+01           \n",
       "17 4.703209e+01      0.000000000             2.467866e+01           \n",
       "18 3.274843e+01      0.000000000             1.444964e+01           \n",
       "19 2.016645e+01      0.000000000             7.017537e+00           \n",
       "20 1.033368e+01      0.000000000             2.880891e+00           \n",
       "21 4.207847e+00      0.000000000             1.007775e+00           \n",
       "22 1.331383e+00      0.000000000             2.932354e-01           \n",
       "23 3.225228e-01      0.000000000             6.916034e-02           \n",
       "24 5.799678e-02      0.000000000             1.290991e-02           \n",
       "25 7.152493e-03      0.000000000             1.863148e-03           \n",
       "26 0.000000e+00      0.000000000             0.000000e+00           \n",
       "27 0.000000e+00      0.000000000             0.000000e+00           \n",
       "28 0.000000e+00      0.000000000             0.000000e+00           \n",
       "29 0.000000e+00      0.000000000             0.000000e+00           \n",
       "30 0.000000e+00      0.000000000             0.000000e+00           \n",
       "   imported_backproj_incl_unobs domestic_backproj_incl_unobs\n",
       "1  0.003141247                  3.056866e+02                \n",
       "2  0.000000000                  2.859795e+02                \n",
       "3  0.000000000                  2.625939e+02                \n",
       "4  0.000000000                  2.398983e+02                \n",
       "5  0.000000000                  2.220498e+02                \n",
       "6  0.000000000                  2.087914e+02                \n",
       "7  0.000000000                  1.973702e+02                \n",
       "8  0.000000000                  1.847542e+02                \n",
       "9  0.000000000                  1.687925e+02                \n",
       "10 0.000000000                  1.508628e+02                \n",
       "11 0.000000000                  1.341985e+02                \n",
       "12 0.000000000                  1.189053e+02                \n",
       "13 0.000000000                  1.028381e+02                \n",
       "14 0.000000000                  8.550350e+01                \n",
       "15 0.000000000                  6.818233e+01                \n",
       "16 0.000000000                  5.200660e+01                \n",
       "17 0.000000000                  3.783090e+01                \n",
       "18 0.000000000                  2.593969e+01                \n",
       "19 0.000000000                  1.598117e+01                \n",
       "20 0.000000000                  8.262887e+00                \n",
       "21 0.000000000                  3.407351e+00                \n",
       "22 0.000000000                  1.092231e+00                \n",
       "23 0.000000000                  2.676521e-01                \n",
       "24 0.000000000                  4.857183e-02                \n",
       "25 0.000000000                  6.030619e-03                \n",
       "26 0.000000000                  0.000000e+00                \n",
       "27 0.000000000                  0.000000e+00                \n",
       "28 0.000000000                  0.000000e+00                \n",
       "29 0.000000000                  0.000000e+00                \n",
       "30 0.000000000                  0.000000e+00                "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now backprojecting unobs + obs by the incubation period\n",
    "\n",
    "# backprojection of imported_unobs cases\n",
    "if (sum(df_cases$imported_unobs)>0) {\n",
    "    sts = new(\"sts\", epoch=df_cases$time, observed=df_cases$imported_backproj_unobs + df_cases$imported)\n",
    "    bpnp.control = list(k = k_used[1], eps = rep(1e-4,2), iter.max=rep(1000,2), Tmark = nrow(df_cases)-10,\n",
    "                        Tmark = tstar - 1, #attention\n",
    "                        B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, \n",
    "                        eq3a.method = c(\"R\",\"C\"))\n",
    "    sts_bp = backprojNP(sts, incu.pmf=inc_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "    # output result\n",
    "    df_cases['imported_backproj_incl_unobs'] = upperbound(sts_bp) %>% as.numeric %>% replace(. < epsilon, 0) \n",
    "    df_cases['imported_backproj_incl_unobs'] = df_cases['imported_backproj_incl_unobs']/sum(df_cases['imported_backproj_incl_unobs'])*sum(df_cases$imported_backproj_unobs + df_cases$imported)\n",
    "} else \n",
    "    df_cases['imported_backproj_incl_unobs'] = 0\n",
    "\n",
    "# backprojection of domestic_unobs cases\n",
    "if (sum(df_cases$domestic_unobs)>0) {\n",
    "    sts = new(\"sts\", epoch=df_cases$time, observed=df_cases$domestic_backproj_unobs + df_cases$domestic)\n",
    "    bpnp.control = list(k = k_used[2], eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                        B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, \n",
    "                        eq3a.method = c(\"C\"))\n",
    "    sts_bp = backprojNP(sts, incu.pmf=inc_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "    df_cases['domestic_backproj_incl_unobs'] = upperbound(sts_bp) %>% as.numeric %>% replace(. < epsilon, 0) \n",
    "    df_cases['domestic_backproj_incl_unobs'] = df_cases['domestic_backproj_incl_unobs']/sum(df_cases['domestic_backproj_incl_unobs'])*sum(df_cases$domestic_backproj_unobs + df_cases$domestic)\n",
    "} else\n",
    "    df_cases['domestic_backproj_incl_unobs'] = 0\n",
    "\n",
    "df_cases %>% filter(time<tstar) %>% tail(30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "342"
      ],
      "text/latex": [
       "342"
      ],
      "text/markdown": [
       "342"
      ],
      "text/plain": [
       "[1] 342"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# total number of imported cases\n",
    "sum(df_cases$imported) + sum(df_cases$imported_unobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "15062"
      ],
      "text/latex": [
       "15062"
      ],
      "text/markdown": [
       "15062"
      ],
      "text/plain": [
       "[1] 15062"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# total number of domestically acquired infections\n",
    "sum(df_cases$domestic) + sum(df_cases$domestic_unobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2889"
      ],
      "text/latex": [
       "2889"
      ],
      "text/markdown": [
       "2889"
      ],
      "text/plain": [
       "[1] 2889"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# among them are unobsered\n",
    "sum(df_cases$domestic_unobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cases %<>% select(-domestic_backproj_unobs,-imported_backproj_unobs) %>% filter(time<tstar)\n",
    "df_cases %>% write.csv(\"../results/JapaneseDataCOVID19_with_backproj (\"%&%format(datestar,\"%y%m%d\")%&%\").csv\", row.names=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color=\"green\">[R] Effective reproduction number</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"purple\">Stan program</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## The code is written by Andrei R. Akhmetzhanov and verified by Kenji Mizumoto\n",
    "\"functions {\n",
    "    /* discretesized version of lognormal distribution */\n",
    "    // the first element is after day one \n",
    "    vector plognormal(real mu, real sigma, int K) {\n",
    "        vector[K] res; \n",
    "        for (k in 1:K)\n",
    "            res[k] = exp(lognormal_lcdf(k | mu, sigma));\n",
    "\n",
    "        return append_row(res[1], res[2:K]-res[1:(K-1)]);        \n",
    "    }\n",
    "\n",
    "    /* discretesized version of Weibull distribution */\n",
    "    vector pweibull(real kappa, real theta, int K) {\n",
    "        vector[K] res; \n",
    "        for (k in 1:K) \n",
    "            res[k] = -expm1(-pow(1.0*k/theta, kappa)); // instead of using Stan function weibull_lcdf it is easier to write it directly\n",
    "\n",
    "        return append_row(res[1], res[2:K]-res[1:(K-1)]);  \n",
    "    }\n",
    "\n",
    "    /* calculating the convolutions */\n",
    "    // X: first function, Yrev: reversed version of the second function\n",
    "    // K: length of X and Yrev\n",
    "    // the result is the vector of length K-1, because the first element equal zero is omitted\n",
    "    vector convolution(vector X, vector Yrev, int K) {\n",
    "        vector[K-1] res;\n",
    "        res[1] = X[1]*Yrev[K];\n",
    "        for (k in 2:K-1) // 2:K-1 is equivalent to 2:(K-1)\n",
    "            res[k] = dot_product(head(X, k), tail(Yrev, k)); // by definition of the convolution\n",
    "\n",
    "        return res;        \n",
    "    }\n",
    "\n",
    "    /* Special form of the convolution with adjusting to the delay */\n",
    "    // F: cumulative distribution function of the reporting delay\n",
    "    // it is the same idea as in [Tsuzuki et al 2017] (doi:10.2807/1560-7917.ES.2017.22.46.17-00710)\n",
    "    vector convolution_with_delay(vector X, vector Yrev, vector F, int K) {\n",
    "        vector[K-1] res;\n",
    "        vector[K] Z = X ./ F; // first we need to nowcast the observed counts (=to account for not yet reported)\n",
    "\n",
    "        res[1] = F[2]*Z[1]*Yrev[K];\n",
    "        for (k in 3:K) \n",
    "            res[k-1] = F[k]*dot_product(head(Z, k-1), tail(Yrev, k-1));\n",
    "\n",
    "        return res;        \n",
    "    }\n",
    "}\n",
    "\n",
    "data {\n",
    "    int<lower = 1> K; //number of days\n",
    "    vector<lower = 0>[K] imported_backproj;\n",
    "    vector<lower = 0>[K] domestic_backproj;\n",
    "    int<lower = K> upper_bound;\n",
    "\n",
    "    // serial interval\n",
    "    real<lower = 0> param1_SI;\n",
    "    real<lower = 0> param2_SI;\n",
    "\n",
    "    // reporting delay is given by Weibul distribution\n",
    "    real<lower = 0> param1_delay;\n",
    "    real<lower = 0> param2_delay;\n",
    "\n",
    "    // incubation period\n",
    "    real mu_inc;\n",
    "    real<lower = 0> sigma_inc;\n",
    "}\n",
    "\n",
    "transformed data {\n",
    "    // RStan requires to declare all variables at the beginning, so there is some mixture\n",
    "    vector[K] cases_backproj;\n",
    "\n",
    "    vector[K-1] conv;\n",
    "    vector[K-1] conv_delay_adj;\n",
    "\n",
    "    cases_backproj = imported_backproj + domestic_backproj;\n",
    "    // all variables declared within {..} are local \n",
    "    // so they will not be recorded in the output of Stan\n",
    "    {\n",
    "        // serial interval\n",
    "        vector[K] gt;\n",
    "        vector[K] gtrev;\n",
    "        // incubation period\n",
    "        vector[upper_bound] IncPeriod; \n",
    "        vector[upper_bound] IncPeriod_inv;\n",
    "        // reporting delay\n",
    "        vector[upper_bound] ReportDelay;\n",
    "        vector[upper_bound] conv_report_inc;\n",
    "        \n",
    "        vector[upper_bound] F_tmp;\n",
    "        vector[K] F; \n",
    "\n",
    "        // serial interval is modeled by Weibull distribution\n",
    "        gt = pweibull(param1_SI, param2_SI, K);\n",
    "        for (k in 1:K) \n",
    "            gtrev[k] = gt[K+1-k];\n",
    "\n",
    "        // incubation period is modeled by lognormal distribution\n",
    "        IncPeriod = plognormal(mu_inc, sigma_inc, upper_bound);\n",
    "        // we need to have it in reverse order for our convolution function used later\n",
    "        for (k in 1:upper_bound)\n",
    "            IncPeriod_inv[k] = IncPeriod[upper_bound+1-k];\n",
    "\n",
    "        // reporting delay is modeled by Weibull distribution\n",
    "        ReportDelay = pweibull(param1_delay, param2_delay, upper_bound);\n",
    "\n",
    "        // convolution of the reporting delay and incubation period\n",
    "        conv_report_inc = append_row(rep_vector(0,1), convolution(ReportDelay, IncPeriod_inv, upper_bound));\n",
    "        conv_report_inc = cumulative_sum(conv_report_inc);\n",
    "        \n",
    "        // here the form is a bit complicated because technically upper_bound is not K in general, it can be any number > K\n",
    "        for (k in 1:upper_bound) \n",
    "            F_tmp[k] = conv_report_inc[upper_bound+1-k];\n",
    "        F = head(F_tmp, K);\n",
    "\n",
    "        // convolution without adjusting for the delay\n",
    "        conv = convolution(cases_backproj, gtrev, K);\n",
    "\n",
    "        // convolution with adjusting for the delay\n",
    "        conv_delay_adj = convolution_with_delay(cases_backproj, gtrev, F, K);\n",
    "    }\n",
    "}\n",
    "\n",
    "parameters {\n",
    "    // effective reproduction number without adjustment to the delay in reporting\n",
    "    vector<lower = 0>[K-1] Rt;\n",
    "\n",
    "    // effective reproduction number with adjustment to the delay in reporting\n",
    "    vector<lower = 0>[K-1] Rt_adj;\n",
    "}\n",
    "\n",
    "model {\n",
    "    Rt ~ normal(2.4, 2.0);\n",
    "    Rt_adj ~ normal(2.4, 2.0);\n",
    "    \n",
    "    // two independent likelihood written in one line\n",
    "    target += gamma_lpdf(tail(domestic_backproj, K-1) | Rt .* conv + 1e-13, 1.0)\n",
    "            + gamma_lpdf(tail(domestic_backproj, K-1) | Rt_adj .* conv_delay_adj + 1e-13, 1.0);\n",
    "}\" %>% cat(file=\"fit_infection.stan\", sep=\"\", fill=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#model discription\n",
    "stanmodel = NULL\n",
    "stanmodel = stan_model(file='fit_infection.stan')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"purple\">Data</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "stan_data = list(\n",
    "  K = nrow(df_cases),\n",
    "  imported_backproj = df_cases$imported_backproj,\n",
    "  domestic_backproj = df_cases$domestic_backproj,\n",
    "  upper_bound = tstar,\n",
    "  ## Serial interval [Nishiura et al 2020 - only certain cases]\n",
    "  param1_SI = 2.305,\n",
    "  param2_SI = 5.452,\n",
    "  ## reported delay\n",
    "  param1_delay = onset2report_fit$param1,\n",
    "  param2_delay = onset2report_fit$param2,\n",
    "  ## incubation period [Linton et al 2020 - with right truncation excl. Wuhan residents]\n",
    "  mu_inc = inc_fit$meanlog,\n",
    "  sigma_inc = inc_fit$sdlog\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"purple\">MCMC settings</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iter =   10000\n",
    "Warmup =  2000\n",
    "Chains =    2\n",
    "Seed =   1234"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"purple\">Model fit</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit <- sampling(stanmodel, \n",
    "                data=stan_data,\n",
    "                iter=Iter,\n",
    "                warmup=Warmup,\n",
    "                chains = Chains, \n",
    "                seed=1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: tibble\n",
      "\n",
      "Loading required package: stringr\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A data.frame: 6 × 12</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>variable</th><th scope=col>time</th><th scope=col>mean</th><th scope=col>se_mean</th><th scope=col>sd</th><th scope=col>2.5%</th><th scope=col>25%</th><th scope=col>50%</th><th scope=col>75%</th><th scope=col>97.5%</th><th scope=col>n_eff</th><th scope=col>Rhat</th></tr>\n",
       "\t<tr><th></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>&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><th scope=row>1</th><td>Rt</td><td>2</td><td>3.803845</td><td>0.01137938</td><td>1.617343</td><td>0.9585255</td><td>2.615784</td><td>3.706837</td><td>4.886268</td><td>7.201542</td><td>20200.72</td><td>0.9999394</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>Rt</td><td>3</td><td>3.802308</td><td>0.01105958</td><td>1.614534</td><td>0.9787390</td><td>2.619588</td><td>3.705229</td><td>4.898468</td><td>7.190752</td><td>21311.65</td><td>0.9999047</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>Rt</td><td>4</td><td>3.827900</td><td>0.01124618</td><td>1.614151</td><td>0.9969482</td><td>2.646266</td><td>3.737979</td><td>4.883448</td><td>7.246082</td><td>20600.51</td><td>0.9999704</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>Rt</td><td>5</td><td>3.831425</td><td>0.01121714</td><td>1.616146</td><td>1.0311345</td><td>2.658927</td><td>3.725188</td><td>4.913039</td><td>7.230205</td><td>20758.55</td><td>0.9999470</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>Rt</td><td>6</td><td>3.847621</td><td>0.01095309</td><td>1.633813</td><td>0.9934183</td><td>2.661992</td><td>3.751999</td><td>4.920421</td><td>7.303961</td><td>22250.04</td><td>0.9999310</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>Rt</td><td>7</td><td>3.833713</td><td>0.01135992</td><td>1.623234</td><td>0.9947558</td><td>2.641717</td><td>3.723736</td><td>4.892645</td><td>7.255641</td><td>20417.93</td><td>0.9999616</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 12\n",
       "\\begin{tabular}{r|llllllllllll}\n",
       "  & variable & time & mean & se\\_mean & sd & 2.5\\% & 25\\% & 50\\% & 75\\% & 97.5\\% & n\\_eff & Rhat\\\\\n",
       "  & <chr> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t1 & Rt & 2 & 3.803845 & 0.01137938 & 1.617343 & 0.9585255 & 2.615784 & 3.706837 & 4.886268 & 7.201542 & 20200.72 & 0.9999394\\\\\n",
       "\t2 & Rt & 3 & 3.802308 & 0.01105958 & 1.614534 & 0.9787390 & 2.619588 & 3.705229 & 4.898468 & 7.190752 & 21311.65 & 0.9999047\\\\\n",
       "\t3 & Rt & 4 & 3.827900 & 0.01124618 & 1.614151 & 0.9969482 & 2.646266 & 3.737979 & 4.883448 & 7.246082 & 20600.51 & 0.9999704\\\\\n",
       "\t4 & Rt & 5 & 3.831425 & 0.01121714 & 1.616146 & 1.0311345 & 2.658927 & 3.725188 & 4.913039 & 7.230205 & 20758.55 & 0.9999470\\\\\n",
       "\t5 & Rt & 6 & 3.847621 & 0.01095309 & 1.633813 & 0.9934183 & 2.661992 & 3.751999 & 4.920421 & 7.303961 & 22250.04 & 0.9999310\\\\\n",
       "\t6 & Rt & 7 & 3.833713 & 0.01135992 & 1.623234 & 0.9947558 & 2.641717 & 3.723736 & 4.892645 & 7.255641 & 20417.93 & 0.9999616\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 12\n",
       "\n",
       "| <!--/--> | variable &lt;chr&gt; | time &lt;dbl&gt; | mean &lt;dbl&gt; | se_mean &lt;dbl&gt; | sd &lt;dbl&gt; | 2.5% &lt;dbl&gt; | 25% &lt;dbl&gt; | 50% &lt;dbl&gt; | 75% &lt;dbl&gt; | 97.5% &lt;dbl&gt; | n_eff &lt;dbl&gt; | Rhat &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| 1 | Rt | 2 | 3.803845 | 0.01137938 | 1.617343 | 0.9585255 | 2.615784 | 3.706837 | 4.886268 | 7.201542 | 20200.72 | 0.9999394 |\n",
       "| 2 | Rt | 3 | 3.802308 | 0.01105958 | 1.614534 | 0.9787390 | 2.619588 | 3.705229 | 4.898468 | 7.190752 | 21311.65 | 0.9999047 |\n",
       "| 3 | Rt | 4 | 3.827900 | 0.01124618 | 1.614151 | 0.9969482 | 2.646266 | 3.737979 | 4.883448 | 7.246082 | 20600.51 | 0.9999704 |\n",
       "| 4 | Rt | 5 | 3.831425 | 0.01121714 | 1.616146 | 1.0311345 | 2.658927 | 3.725188 | 4.913039 | 7.230205 | 20758.55 | 0.9999470 |\n",
       "| 5 | Rt | 6 | 3.847621 | 0.01095309 | 1.633813 | 0.9934183 | 2.661992 | 3.751999 | 4.920421 | 7.303961 | 22250.04 | 0.9999310 |\n",
       "| 6 | Rt | 7 | 3.833713 | 0.01135992 | 1.623234 | 0.9947558 | 2.641717 | 3.723736 | 4.892645 | 7.255641 | 20417.93 | 0.9999616 |\n",
       "\n"
      ],
      "text/plain": [
       "  variable time mean     se_mean    sd       2.5%      25%      50%     \n",
       "1 Rt       2    3.803845 0.01137938 1.617343 0.9585255 2.615784 3.706837\n",
       "2 Rt       3    3.802308 0.01105958 1.614534 0.9787390 2.619588 3.705229\n",
       "3 Rt       4    3.827900 0.01124618 1.614151 0.9969482 2.646266 3.737979\n",
       "4 Rt       5    3.831425 0.01121714 1.616146 1.0311345 2.658927 3.725188\n",
       "5 Rt       6    3.847621 0.01095309 1.633813 0.9934183 2.661992 3.751999\n",
       "6 Rt       7    3.833713 0.01135992 1.623234 0.9947558 2.641717 3.723736\n",
       "  75%      97.5%    n_eff    Rhat     \n",
       "1 4.886268 7.201542 20200.72 0.9999394\n",
       "2 4.898468 7.190752 21311.65 0.9999047\n",
       "3 4.883448 7.246082 20600.51 0.9999704\n",
       "4 4.913039 7.230205 20758.55 0.9999470\n",
       "5 4.920421 7.303961 22250.04 0.9999310\n",
       "6 4.892645 7.255641 20417.93 0.9999616"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "require(tibble)\n",
    "require(stringr)\n",
    "\n",
    "# Check the result\n",
    "## Rt is for backproj\n",
    "## Rt_adj is for backproj_incl_unobs\n",
    "## The first and the last values should be ignored because no information is available and Rt* is unstable\n",
    "### notice that Reff starts from t=2 not t=1, so one of the coordinates is shifted by one\n",
    "summary(fit)$summary %T>% \n",
    "    write.csv(\"../results/Results of Stan (script C).csv\") %>%\n",
    "    as.data.frame %>% tibble::rownames_to_column(\"index\") -> df_summary\n",
    "\n",
    "out = stringr::str_match_all(df_summary$index, \"[\\\\d+(.*)]\")\n",
    "df_summary$time = sapply(out, function(x) paste(x[,1],collapse=\"\"))\n",
    "df_summary %<>% mutate(time = as.numeric(time) + 1, variable = str_extract(index, \"^[^\\\\[]+\")) %>%\n",
    "    select(-index) %>% select(variable, time, everything())\n",
    "\n",
    "df_summary %>% head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A data.frame: 6 × 6</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>time</th><th scope=col>Rt_lower</th><th scope=col>Rt_IQR_lower</th><th scope=col>Rt_median</th><th scope=col>Rt_IQR_upper</th><th scope=col>Rt_upper</th></tr>\n",
       "\t<tr><th></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><th scope=row>1</th><td>2</td><td>0.9606674</td><td>2.661039</td><td>3.735816</td><td>4.905539</td><td>7.243792</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>3</td><td>0.9695339</td><td>2.637299</td><td>3.716757</td><td>4.891109</td><td>7.245668</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>4</td><td>0.9545927</td><td>2.646601</td><td>3.737445</td><td>4.894425</td><td>7.228038</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>5</td><td>0.9811007</td><td>2.649600</td><td>3.756910</td><td>4.899631</td><td>7.293788</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>6</td><td>1.0486528</td><td>2.667259</td><td>3.752073</td><td>4.909372</td><td>7.228146</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>7</td><td>1.0087498</td><td>2.690629</td><td>3.757193</td><td>4.898231</td><td>7.266286</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 6\n",
       "\\begin{tabular}{r|llllll}\n",
       "  & time & Rt\\_lower & Rt\\_IQR\\_lower & Rt\\_median & Rt\\_IQR\\_upper & Rt\\_upper\\\\\n",
       "  & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t1 & 2 & 0.9606674 & 2.661039 & 3.735816 & 4.905539 & 7.243792\\\\\n",
       "\t2 & 3 & 0.9695339 & 2.637299 & 3.716757 & 4.891109 & 7.245668\\\\\n",
       "\t3 & 4 & 0.9545927 & 2.646601 & 3.737445 & 4.894425 & 7.228038\\\\\n",
       "\t4 & 5 & 0.9811007 & 2.649600 & 3.756910 & 4.899631 & 7.293788\\\\\n",
       "\t5 & 6 & 1.0486528 & 2.667259 & 3.752073 & 4.909372 & 7.228146\\\\\n",
       "\t6 & 7 & 1.0087498 & 2.690629 & 3.757193 & 4.898231 & 7.266286\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 6\n",
       "\n",
       "| <!--/--> | time &lt;dbl&gt; | Rt_lower &lt;dbl&gt; | Rt_IQR_lower &lt;dbl&gt; | Rt_median &lt;dbl&gt; | Rt_IQR_upper &lt;dbl&gt; | Rt_upper &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|\n",
       "| 1 | 2 | 0.9606674 | 2.661039 | 3.735816 | 4.905539 | 7.243792 |\n",
       "| 2 | 3 | 0.9695339 | 2.637299 | 3.716757 | 4.891109 | 7.245668 |\n",
       "| 3 | 4 | 0.9545927 | 2.646601 | 3.737445 | 4.894425 | 7.228038 |\n",
       "| 4 | 5 | 0.9811007 | 2.649600 | 3.756910 | 4.899631 | 7.293788 |\n",
       "| 5 | 6 | 1.0486528 | 2.667259 | 3.752073 | 4.909372 | 7.228146 |\n",
       "| 6 | 7 | 1.0087498 | 2.690629 | 3.757193 | 4.898231 | 7.266286 |\n",
       "\n"
      ],
      "text/plain": [
       "  time Rt_lower  Rt_IQR_lower Rt_median Rt_IQR_upper Rt_upper\n",
       "1 2    0.9606674 2.661039     3.735816  4.905539     7.243792\n",
       "2 3    0.9695339 2.637299     3.716757  4.891109     7.245668\n",
       "3 4    0.9545927 2.646601     3.737445  4.894425     7.228038\n",
       "4 5    0.9811007 2.649600     3.756910  4.899631     7.293788\n",
       "5 6    1.0486528 2.667259     3.752073  4.909372     7.228146\n",
       "6 7    1.0087498 2.690629     3.757193  4.898231     7.266286"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_summary %>%\n",
    "    filter(variable=='Rt_adj') %>% \n",
    "    select(time,`2.5%`,`25%`,`50%`,`75%`,`97.5%`) -> df_summary_\n",
    "names(df_summary_) = c('time','Rt_lower','Rt_IQR_lower','Rt_median','Rt_IQR_upper','Rt_upper')\n",
    "df_summary_ %>% head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 6 × 15</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>date</th><th scope=col>time</th><th scope=col>domestic_unobs</th><th scope=col>imported_unobs</th><th scope=col>domestic</th><th scope=col>imported</th><th scope=col>imported_backproj</th><th scope=col>domestic_backproj</th><th scope=col>imported_backproj_incl_unobs</th><th scope=col>domestic_backproj_incl_unobs</th><th scope=col>Rt_lower</th><th scope=col>Rt_IQR_lower</th><th scope=col>Rt_median</th><th scope=col>Rt_IQR_upper</th><th scope=col>Rt_upper</th></tr>\n",
       "\t<tr><th scope=col>&lt;date&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><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>2020-02-02</td><td>39</td><td>0</td><td>0</td><td>4</td><td>0</td><td>1.1739613</td><td>2.409345</td><td>1.1749413</td><td>2.728569</td><td>0.2426307</td><td>0.6056138</td><td>0.8556800</td><td>1.1559016</td><td>1.768350</td></tr>\n",
       "\t<tr><td>2020-02-03</td><td>40</td><td>0</td><td>0</td><td>9</td><td>1</td><td>0.8989732</td><td>1.507233</td><td>0.9008537</td><td>2.270909</td><td>0.1165546</td><td>0.3413098</td><td>0.5030566</td><td>0.6928784</td><td>1.129236</td></tr>\n",
       "\t<tr><td>2020-02-04</td><td>41</td><td>0</td><td>0</td><td>5</td><td>0</td><td>0.5000108</td><td>1.472679</td><td>0.5019078</td><td>3.032050</td><td>0.1035112</td><td>0.3041492</td><td>0.4552194</td><td>0.6290248</td><td>1.022204</td></tr>\n",
       "\t<tr><td>2020-02-05</td><td>42</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0.2010538</td><td>2.091894</td><td>0.2022331</td><td>4.879109</td><td>0.1479826</td><td>0.4044045</td><td>0.5753616</td><td>0.7828775</td><td>1.226896</td></tr>\n",
       "\t<tr><td>2020-02-06</td><td>43</td><td>0</td><td>0</td><td>5</td><td>1</td><td>0.1082749</td><td>3.269345</td><td>0.1083783</td><td>7.486603</td><td>0.2899271</td><td>0.6441551</td><td>0.8832466</td><td>1.1554194</td><td>1.716245</td></tr>\n",
       "\t<tr><td>2020-02-07</td><td>44</td><td>0</td><td>0</td><td>1</td><td>2</td><td>0.1967178</td><td>4.742397</td><td>0.1943830</td><td>9.907324</td><td>0.5015735</td><td>1.0310168</td><td>1.3605620</td><td>1.7081972</td><td>2.470284</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 6 × 15\n",
       "\\begin{tabular}{lllllllllllllll}\n",
       " date & time & domestic\\_unobs & imported\\_unobs & domestic & imported & imported\\_backproj & domestic\\_backproj & imported\\_backproj\\_incl\\_unobs & domestic\\_backproj\\_incl\\_unobs & Rt\\_lower & Rt\\_IQR\\_lower & Rt\\_median & Rt\\_IQR\\_upper & Rt\\_upper\\\\\n",
       " <date> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\t 2020-02-02 & 39 & 0 & 0 & 4 & 0 & 1.1739613 & 2.409345 & 1.1749413 & 2.728569 & 0.2426307 & 0.6056138 & 0.8556800 & 1.1559016 & 1.768350\\\\\n",
       "\t 2020-02-03 & 40 & 0 & 0 & 9 & 1 & 0.8989732 & 1.507233 & 0.9008537 & 2.270909 & 0.1165546 & 0.3413098 & 0.5030566 & 0.6928784 & 1.129236\\\\\n",
       "\t 2020-02-04 & 41 & 0 & 0 & 5 & 0 & 0.5000108 & 1.472679 & 0.5019078 & 3.032050 & 0.1035112 & 0.3041492 & 0.4552194 & 0.6290248 & 1.022204\\\\\n",
       "\t 2020-02-05 & 42 & 0 & 0 & 1 & 1 & 0.2010538 & 2.091894 & 0.2022331 & 4.879109 & 0.1479826 & 0.4044045 & 0.5753616 & 0.7828775 & 1.226896\\\\\n",
       "\t 2020-02-06 & 43 & 0 & 0 & 5 & 1 & 0.1082749 & 3.269345 & 0.1083783 & 7.486603 & 0.2899271 & 0.6441551 & 0.8832466 & 1.1554194 & 1.716245\\\\\n",
       "\t 2020-02-07 & 44 & 0 & 0 & 1 & 2 & 0.1967178 & 4.742397 & 0.1943830 & 9.907324 & 0.5015735 & 1.0310168 & 1.3605620 & 1.7081972 & 2.470284\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 6 × 15\n",
       "\n",
       "| date &lt;date&gt; | time &lt;dbl&gt; | domestic_unobs &lt;dbl&gt; | imported_unobs &lt;dbl&gt; | domestic &lt;dbl&gt; | imported &lt;dbl&gt; | imported_backproj &lt;dbl&gt; | domestic_backproj &lt;dbl&gt; | imported_backproj_incl_unobs &lt;dbl&gt; | domestic_backproj_incl_unobs &lt;dbl&gt; | Rt_lower &lt;dbl&gt; | Rt_IQR_lower &lt;dbl&gt; | Rt_median &lt;dbl&gt; | Rt_IQR_upper &lt;dbl&gt; | Rt_upper &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| 2020-02-02 | 39 | 0 | 0 | 4 | 0 | 1.1739613 | 2.409345 | 1.1749413 | 2.728569 | 0.2426307 | 0.6056138 | 0.8556800 | 1.1559016 | 1.768350 |\n",
       "| 2020-02-03 | 40 | 0 | 0 | 9 | 1 | 0.8989732 | 1.507233 | 0.9008537 | 2.270909 | 0.1165546 | 0.3413098 | 0.5030566 | 0.6928784 | 1.129236 |\n",
       "| 2020-02-04 | 41 | 0 | 0 | 5 | 0 | 0.5000108 | 1.472679 | 0.5019078 | 3.032050 | 0.1035112 | 0.3041492 | 0.4552194 | 0.6290248 | 1.022204 |\n",
       "| 2020-02-05 | 42 | 0 | 0 | 1 | 1 | 0.2010538 | 2.091894 | 0.2022331 | 4.879109 | 0.1479826 | 0.4044045 | 0.5753616 | 0.7828775 | 1.226896 |\n",
       "| 2020-02-06 | 43 | 0 | 0 | 5 | 1 | 0.1082749 | 3.269345 | 0.1083783 | 7.486603 | 0.2899271 | 0.6441551 | 0.8832466 | 1.1554194 | 1.716245 |\n",
       "| 2020-02-07 | 44 | 0 | 0 | 1 | 2 | 0.1967178 | 4.742397 | 0.1943830 | 9.907324 | 0.5015735 | 1.0310168 | 1.3605620 | 1.7081972 | 2.470284 |\n",
       "\n"
      ],
      "text/plain": [
       "  date       time domestic_unobs imported_unobs domestic imported\n",
       "1 2020-02-02 39   0              0              4        0       \n",
       "2 2020-02-03 40   0              0              9        1       \n",
       "3 2020-02-04 41   0              0              5        0       \n",
       "4 2020-02-05 42   0              0              1        1       \n",
       "5 2020-02-06 43   0              0              5        1       \n",
       "6 2020-02-07 44   0              0              1        2       \n",
       "  imported_backproj domestic_backproj imported_backproj_incl_unobs\n",
       "1 1.1739613         2.409345          1.1749413                   \n",
       "2 0.8989732         1.507233          0.9008537                   \n",
       "3 0.5000108         1.472679          0.5019078                   \n",
       "4 0.2010538         2.091894          0.2022331                   \n",
       "5 0.1082749         3.269345          0.1083783                   \n",
       "6 0.1967178         4.742397          0.1943830                   \n",
       "  domestic_backproj_incl_unobs Rt_lower  Rt_IQR_lower Rt_median Rt_IQR_upper\n",
       "1 2.728569                     0.2426307 0.6056138    0.8556800 1.1559016   \n",
       "2 2.270909                     0.1165546 0.3413098    0.5030566 0.6928784   \n",
       "3 3.032050                     0.1035112 0.3041492    0.4552194 0.6290248   \n",
       "4 4.879109                     0.1479826 0.4044045    0.5753616 0.7828775   \n",
       "5 7.486603                     0.2899271 0.6441551    0.8832466 1.1554194   \n",
       "6 9.907324                     0.5015735 1.0310168    1.3605620 1.7081972   \n",
       "  Rt_upper\n",
       "1 1.768350\n",
       "2 1.129236\n",
       "3 1.022204\n",
       "4 1.226896\n",
       "5 1.716245\n",
       "6 2.470284"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_cases %>%\n",
    "    mutate(date = datemin + time) %>% \n",
    "    left_join(df_summary_, by='time') %>%\n",
    "    select(date,time,everything()) %>% na.omit -> figure_data\n",
    "\n",
    "range = c(as.Date(\"2020-02-01\"), datestar - 17)\n",
    "figure_data %<>% filter(date < range[2], date > range[1])\n",
    "\n",
    "figure_data %>% head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAAJYCAMAAABvmDbGAAAAUVBMVEUAAAATgKEzMzNNTU1Z\nWVlkqLxoaGhpkG5rr8R8fHyMjIyampqnp6eysrK10Ni9vb3AoDrE3+fHx8fQ0NDZ2dnh4eHp\n6enr6+vw8PD6qxj///+L0ViQAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2diXaj\nvBKE8e8bB8dOYrJ4HN7/Qa/ZtbSkBiSQoOqcGdtQaIH+IiFAZCUEQbOVrV0ACNqCABIEeRBA\ngiAPAkgQ5EEACYI8CCBBkAcBJAjyIIAEQR4EkCDIgwASBHlQGJAyQdKKx6VZzUkiSMksSVeL\np+YasLRQEloYpOYnQIK2plAgjV0xwzlaFpA8JwntRgDJS44Aae9aAqTnr+sxO323Xb52dZY9\nTtnb88vPOcvO34Ykrqfnyp92weexN0rL29R79Slesvqc7Fp/fD+XXh6l1IdrPp5rju+lUDIh\nxWrdp1AhYUPBVy+9H0/v8sa/l2N2vPw+S9QU5Ceri3zJfohS9+UTa/dMInsm8SPXS14MRaBl\nQDrXBH0rIL1l2bUsv5pTqSuZxDFrt6wWvA9Gafm5/95ISPGYPaP4nh3LiqZKx1IHqVlzEUo2\npNjmaQCp99VLn5jLG39nfTGP9WbXpvTPUuilHson1u5nSEKsl7gYikHLgHT8KR9vTZj1q5+h\nVP35/a3/9P6e1ZioPe/Ze/3/uUnl6xmbNRrK8i71RmKKP1Wb91a1A7/Z8bt8nKswVED6zbIm\n4QGkPkVxnVS1xizV66wW57nx+6N8PAm5P8ta1e9Yg/KdvVOl7ssn1u5U5f9E6CTXS1gMRaHg\no3b1ryqIHkP8tkvrjsk1q3szj7qXJyVRVvHyEPz3esXxqi7vU28kpXjJvptu1SX7rBeeNJCu\n9Zryccz0FNt1XwaQxHpds+Nd3bhpZi/Pz3tViHt2qipxef6nlVoon1y7fo+I9cI5WWxaBqR2\noQJSvfREj5P3ofL7/X5u/O0f77eTsjyT/GqKx6eapQ8xaWGrbs1bJpdMWPcwgCR8/azbLmXj\nBv17Rcc5qxqZr6qxGZpmqdR9+cTavT3PhL7uWr2ExVAUWmiwofvQQTJccGpXfx4HHNtuzClT\nlqshKaf407Z8QvoaBM3iMwFSt84N0lMnen3z5evZBD2L/jR9V8QRpR7qLtTuXn89fSr1EhZD\nUSgCkCxJPP/Mn69fdwmkKt6l5baQHAOScI5ErjNUpjtd+urGPrSN21QvVafsLXtcqrbHVmqx\nds8zqkvFzLtar34xFIVWB6kePJD0NnhOAgTdSdTbSVmuhqScYte1Oxq7du3X4RxuUtfu/jzL\nOakbC1276nTreRr0/FefshGl7st3UtH/vVSDFNqeahZDUWh1kNpLPb/D+FU9qvAQovK79Tcx\nfbyqy5UspRSf5/Df7WBDc3Z/7Ary6FJ4awL0kwKpHWz4FEHqN1Tq9ajbPmnjS1ugZsT6/Mzo\nNzupDVdXaql8Xe2EXajvKWNzDi2uxUG6l3K4/tYXRn6Pw/B3HTBN8J2qOG4Hn7NqfLj68VCX\nK1mKKf5kXfv1nR1/heHvc/b2aFP4rMfVvzIKpGb4+0s4gxM2VP9AXIXxtPpE5lmQZvi7JvXU\nVaOkSi2UT6xdM85dpyzWS1gMRaFQIIlnxlJvRw3X7qrlcEH2t7kyWRH32SVS/6l/669Bysv7\n1EstxTrqPoULsufW2lzRbC60NpdG3ymQ2guyAkjChipIj6rzJW4sXJAtuwtDb+01Wa3UQ/nE\n2rVXXuu9IdRLXAzFoIVB+j31fZc+hu7VvTLi+FN1J8ylCZHP6habn+/2T/0z1i6/xPKylEJy\nSPFSB+azz1V1ij6feQ9B/HPqbgt6tjjaLUJDitWNO28/QurDhlqX9fNZHmnj/hahsi7Ge+25\nl2oejfryibUrf+p7ge7qnhIXQxEonU72mucD2uViCJIFkByZVucizyYSV2wgqwCSVd0p0tlt\nhXYtgGTXdzXCcf5aI2soJaUDEgRFLIAEQR4EkCDIgwASBHkQQIIgDwJIEORBAAmCPAggQZAH\nASQI8qAQIP0HQUnJQ9AHAUn5/Y+9Jd8ZJFE49+kESHDC6cEJkOCE04MTIMEJpwcnQIITTg9O\ngAQnnB6cAAlOOD04ARKccHpwAiQ44fTgBEhwwunBuRpI/cT0yqehXAAJzqida4F0bP9TP03l\nAkhwRu0ESHDC6cG5KkglQIJzI871QGrOiQiQmrvS/0FQQloNpGOJrh2c23HiHAlOOD04ARKc\ncHpwAiQ44fTgBEhwwunBiTsb4ITTgxP32sEJpwcnQIITTg9OgAQnnB6cAAlOOD04ARKccHpw\nAiQ44fTgBEhwwunBCZDghNODEyDBCacHJ0CCE04PToAEJ5wenAAJTjg9OAESnHB6cAIkOOH0\n4ARIcMLpwQmQ4ITTgxMgwQmnBydAghNOD06ABCecHpwACU44PTgBEpxwenACJDjh9OAESHDC\n6cEJkOCE04MTIMEJpwcnQIITTg9OgAQnnB6cAAnOpZx/tdbKPbATIMG5lBMgOQSQ4OQ4AZJD\nAAlOjhMgOQSQ4OQ4AZJDAAlOjhMgOQSQ4OQ4AZJDAAlOi7PnByA5BJDgtDgBElcACU6LEyBx\nBZDgtDgBElcACU6LEyBxBZDgtDgBElcACU6LEyBxBZDgtDgBElcACU6LUwNJACp87ks5ARKc\ngZ0AiSuABKfFCZC4AkhwWpwAiSuABKfFCZC4AkhwWpwAiSuABKfFCZC4AkhwWpwAiSuABKfF\nCZC4AkhwWpwAiSuABKfFCZC4AkhwWpwAiSuABKfFCZC4AkhwWpwAiSuABKfFCZC4AkhwWpwA\niSuABKfFCZC4AkhwWpwAiSuABCfh1LgBSA4BJDgJJ0AaK4AEJ+EESGMFkOAknABprAASnIQT\nII1V5CAVBdvKTxROpxMgjVX0IJlISuUgpekESGMVP0gGklI5SGk6AdJYJQASjVIqBylNJ0Aa\nqyRAokhK5SCl6QRIYxU3SEVhJCmVg5Sm0w0SB6iYamTTjkDSSUrlIKXpBEhjFQSkf770Mchb\nmhBDDS/Dp2HB2sX0pGhBUn77aJG0JimVv3ZpOtEijRVAgpNwAqSxihqkAiCt5ARIY5UMSCpJ\nqRykNJ0AaawAEpyEEyCNFUCCk3ACpLECSHASToA0VjGDVACktZwAaazSAamwWfmJwslxAqSx\nAkhwEk6ANFYACU7CCZDGCiDBSTgB0lhFDFIBkFZzAqSxSgikwmzlJwonywmQxgogwUk4AdJY\nASQ4CSdAGiuABCfhBEhjBZDgJJwAaaziBUnjSCYplYOUphMgjRVAgpNwAqSxAkhwEk6ANFYA\nCU7CCZDGKimQCoOVnyicPCdAGqtoQaI4AkhLOQHSWAEkOAWnnRuAZBZAglNwAqSpAkhwCk6A\nNFVpgVSQVn6icDqcAGmqYgWJ5gggBXYCJJsyUeo6dm58AaRkneNB+rMQFUONOAJIIxOF0+UE\nSDY19NT/AyQ4bU6AZNOGQCoIKz9ROF1OgGRTgiCZOAJIYZ0AySaABCfTCZBsAkhwMp0AySaA\nBCfTCZBs2hJIhWblJwqn0wmQyvJo9An0ACQ4bU6AVB7NINkEkOAUnADpaGmROmnNUQmQ4JSc\nuwfpaO/a9Tre1XXs3PiaD5KZo56kVA5SWs6dglSp+coFKTur69i58QWQknXuFKT+27G0gdTp\n2bVLZLABIK3k3DlIx/4/qwASnA7n3kFq5PJXIJ3UZezc+AJIyTp3DlKtvQw2FFMThdPtBEj7\nGWwASADJs5MLUqdtnCMBJIDk2Tl28pNkQLJyBJAAkmfnhFmE0hhsAEhrOQGSVY/rKctO79Sq\nBEEqpiUKJ8MJkGy6H9shu4e+DiDBKTgBkk2X7FqdHF2zi74OIMEpOAGSTdUIQ/dPW8fOjS+A\nlKwTINkEkOBkOgGSTefsa1Ndu2JSonAynADJpp/6AlKWnRIZbHBxBJCCOQGSVb/nJ0hvn9Qq\ngASn4ARIUwWQ4BSc00GigIqhRhwBpBGJwslxAiSr7vWdDVfiFClNkIoJicLJcQIkm36PhoeR\nSoAEp+QESDads6oxelyyN30dQIJTcAIkm7oLsalckAVIqzkBkk1Z9lt9/GTEs38RguTmCCCF\ncgIkmy7dLPpXfV2aIBVr7/qtOgGSTfdTAxJxigSQ4BSdAGmqABKcghMgcZTIYANAWs8JkGz6\nFKbjUm+4A0hwcvAASKU8r506cJcoSEUyBykJJ0BiqZvMOI2uHYsjgOTVCZBmaypIzbz97YTj\nysTjACk1J0CarYkgHfsXyRz1l2EApNScAIml4QzJ10yrxxIgbckJkFjyDtKxXB2kIpWDlIQT\nII2TrzkbzCDVb+P8798cffA0Kw9IVkNB90VY8GdfQK1IUCNB+n7zNGp3LNEibcuJFomtx2f1\nkKy+fAJIPTcAaStOgMTUzyXLjlfiAdlJIHUv2gwCEpMjgOTTCZBYqhqjyw+9bs51pFVB+hiT\nKJx2ASSWsiyjGqNmHTs3WQBpS06fIAlEJVH3kg/S17NFugZokYLc2QCQVnACJKZ+3jyeIzkF\nkFJzAiS2Hu9HX6N2Ti0DUjEiUTjtAkhj9H1K4e5vgLSCEyCN06++CCDBCZC4xoRexqzgklcC\nSKGdAImllF7GTIFEwvRRcElK5XACpLWcW3wZswkkjSSA5NEJkFhK6R2yAGkNJ0BiKV2QcoC0\niBMgsZTSy5iNIKkkPUFikpTK4QRIazm3+DJmgLSGEyDxlNDLmAHSGk6ANFuRgVSYQcoBUjAn\nQJqtqEHKAdIyToA0W0mDxCMplcO5MZAEoFao0SgnQOIlys9+r06ANFsJgZQDpFBOgDRbAAlO\ngMR2mhUzSLkTJBZJSYTyuk6AxJIwZbG+jp0bX9NB+rCClAOkQE6AxBJAYmq3ToA0QvfrOfon\nZAHSOk6ANEraey/LqEHSOJJIaq2MRJMI5XWdAGmULmd9GUCCEyCxnWYBJDgBEteYzOQnAGkd\nJ0BiadTkJ59vWVZSwxJ8hQMp10BikJREKK/rBEgsjZj85HFq35BpmCp8UrmmgURxBJDCOAES\nSyPmbGiZ+8qIcYnJ5WLXqABI6zgBEksjQLJZJ5fLI0g5QArhBEgsjZj8BCDt0hkUJA5RaYA0\nYvIT2+nU5HIFBclNUhKhvK4TIPHEn/zk0Q3wGV/xN6Fck0CiOQJIQZwAaba0Ltx7dcnpSrRd\nfIUEKQdIAMm3c3N3NkwCyUlSEqG8rhMgOZVJItazc+MLIKXmBEhOjQXprV6QnWI9R+pJAkge\nnQCJpfv5fK9GHKiROAWkawNbtsqonUiHkSOA5NE5Fw++cxMgnRtaHtmbvk4B6djcG/S7ynWk\niSC5SIo6lNd1AqRa/DsbHu0nsU61lkYrW2FBygESQPLr5N+0WnXt7mfqBjoFmLfs8qgeu1jl\nXjuAtLQTINUa/RgFMYKQGaxznqPwAJKFIx0kB0lRh/K6ToBUi38d6f1susqqduGaZwCvcwbt\nAFIyToBUa2sXZJkg5QDJlxMg1QJIAGmeEyDVCjFB5PVotE4uV3CQ7CRFHcrrOgFSrQAgXS3W\nyeVi1qjggpQDJE/ORUGyAZUGSI1YM60eM/Jhi3ECSKk4AVKtADOtzmqJOs0HSSXncDhYQbKS\nFHUor+sESLUCzLT6ls16EqmRX5AOrQCSfydAqhVg1O5+PM+6hFTLI0gHQTJJAMmDEyDVGjnT\nKueCrG1cYnK5poM0MASQQjgBUq0AtwhFCBJ1lqSBZCMp6lBe1wmQao2+adU9HZcXeQNJ6NBp\nfTuA5MEJkGqNmSBS/JTWsXPjazZI1GAdQArgBEi1goC02iT6xSyQLCRFHcrrOgFSrQBduxUn\n0XeDdABInp0AqVaAwYYVJ9F3gKQ2SQDJgxMg1Qoy/L3a3N86SAcVpANA8utcBSQKqERAsihy\nkIQ2KLeDZCYp6lBe1wmQagUAacVJ9AeQ6M6c8hMgeXACpFr8rp15Qu94JtHngHSwgGQkKepQ\nXtcJkGoFGGxYcRJ9DSTlGqwCFkDy4ARItbggvfGn4/IifyDlKkgHC0gmkqIO5XWdAKlWgAki\nvWgSSIUbJGnJC0Ca7wRItfgtUv3hnLLYMd/+5HItBZKBpKhDeV0nQKrFHmx4a7p2BEfJgXQA\nSADJt3P85Cc6INprXczMTS7XSJBMYw0yWy85l6SoQ3ldJ0CqFQCk7lFzqhc4uVyTQVI5Akie\nneuCJBCVBkg20ZOfPNbs2tlBOgAkgOTZGQCkczdSHmeLlIsgcUmKOpTXdQKkWmyQLNOnmt5G\nMevOhn8T9NHrpVbV+LxokpZ+6JqS847VhPLw+U/61FbMcP5ZVqwuLki26VPpt1G8L39ng9oi\nkQ2SsJhukagmKeo2YV0nWqRaXJCO2ZdxXSwXZIseJFvPTgaJ2beLOpTXdQKkWmMfNSfXsXPj\nKzBIBxtIBElRh/K6ToBUiwvSu2X6VPmC7HrTcXFBygESQPLtZA82vJunT40VJOpyLAsknaSo\nQ3ldJ0CqFeL9SD7kCSSSo56wCiRekxR1KK/rBEi1dglSDpAAkmdniFdfNneKZ6eln5DtQXKc\nIrlB0kiKOpTXdcYBkgDUgnUXFACka9NqZYvP2aCAZDxF6lfVILFIijqU13UCpFpskO7s6biO\nzcyQv+sNNrgapG6dBSSFpKhDeV0nQKrFBel3zNso5M9JmgBS4RUkmaSoQ3kVpzc8+M5NgHTO\nqsbocXE9IVtWp0iXynpdeqbVkSAdOpA4TVKMobyuEyDJCjCJfn/T6pxZ9AODlLtBkkiKMZTX\ndQIkWXyQaix+GC9j7mY3nvX+y/kgWcYaeCCJJMUYyus6AZIs/tso2qG4q74ukutIPUiMBqnB\nrAWJQVKMobyuEyDJYk8QeWpAck1+4kuhQco5IA0kxRjK6zoBkqwQF2QtzwCyBZAidwIkWWNB\nYgw22J4BZGsJkA5OkHqUYgzldZ0ASRYXpE9hEqFPZZ12QVY1TNBskBwc1YYOJBtJBZF9scH3\nVoyd9xwgyZoyHZc6cGe4IDtL40HqAt83SHWE/dPyMRUjLjy4TluVABJH7EfNW3oYXbs3yzOA\nbM0AydCzK4Z3vfR9O3GlFaV/SiYbA8lap7RA4hAVwzkSJe2CrPkZQLZ8g1Sv0pokLkj6u/22\nNN2+o04AiaMgE0SuMthgA6lbB5AIp6tSAImjCQ/2uaYsjg0kIU4m9+0IkDYy3b67UgCJowAg\nedFckARIxDiZ2iTtBiT2xC9EKN9uN4A0RwmDxCWJAmkbs4S7a8UF6VapAEjTFRdILzaQpjZJ\nJEibmCXcF0jFrQGpKACSVQ/eE7Li9aZ1zpEA0hgno1YskIoBpGcHDyAZxX2r+Wog9XGggqSE\nybS+HQ3SBiY3ZtSKA9JNAukGkIy6dO9qIWY0iaJrNwWkF4DEqJYbpBogEaQbQDJpxBOyXjQT\nJFPPTibphbzYNAKk5Cc35lTLCdJNB+lWAiRamwWJ1yTtDCTGY/ZD5N4AUrnZrp2xZyeR9MIe\nbjCBlPqcrMYKO9PsA/ZGglQAJFrcwQZfWgEkG0m7A8n5mH0bsANACkgFQKLFHP72pckg5U6Q\nChkkVpNkBCntyY0tNXY9HdwEbAGQWm3mgqwEkuUUSQaJ2yTtEKRCcmqq47XnhgCpAEhjFQNI\nAyUySGSIyCBxmiSApOqvHfU2g3RbFSQbUJHctKqvY+fG11IgMYcbzCClPN2+laO+bgaQJG5q\nkF6fAkgObQqkwgCSmaSdglRXj0pTbYAqkCqOXp87UbpVCCCRup9jHbWbAxKDpP2CZKq7ClL+\n2oH0+ircKgSQaD0Yk+h7kSeQTMExgMRrkgCSKgWkXATpFSA5FXvXThmKMwWHCSTTBhaQ0n2T\nkrpHZoD0WvHz3JMtSK/9PXcAidYnZxJ9HwoMUiGB5CZp8yCx604N0tUtUcVRvS9fhUEHgKSq\nH2uIdBL96SCxmqStg2SrvBOkpkt3GFT9bh0ASVE3NyTBUQwgCRHBA6kwgcQJpm2ApGFkqLwL\npObcqAZoIOmZFEAaq3hBsoS/BJKTJBtIqb7bj+CIrLwDpGaQoW2IXtuGqQfpBpD4ShwkTpO0\nZZDy3F57O0i5xFE92ND+ZQJIYxUbSJyeXR0/A0gukrYLUq5qHEjNIF3PUTNq1+zQ18YBkNhK\nHSRtnvBdgDTsBgdKNpB6bg7SdaR2AUAap4hAyiWQbNHf9gKZTZIVpDRf22wGSam/BaS+J3do\nb2hoyOqbKIA0SsmD5G6Stg+S+RqAEaSmAepG6V6l0YcBpMIR9MYVAMmHFgCpkECyN0nbA+lD\n5ch4q5QBpNfh8lH9Rb6wlHckASRZcd/9PbAhTQ3EBsnZJG0dpO6uBGoHkCB1p0T9sLdyq0O7\nokYOIAlKESQHR90UeFSTZD/htpGUJEh19U0kESC9dveo9sN1ys13fVO1PkgEUOt37WJ9jGI2\nSOQL/nYAkvRnxEBSV/dmcTc4J3OkPSHbJtesAEiaIn2MwgtINpLkYNLSZheUX6XQTgGkvu40\nSR9ivSWQumFv7S7W5gJTTRlAohV31y4fB9KHwI0GEnXCbeKMXVB+lQI7iwEkAR/9z0kujG2q\nIHXD3iRIRftcxfBc0m5AOh6PxDMSsiJ9jIIEycmRCpKFJLEjZCcpOZCkih+I3WACqbW21131\nyU9ax+5AOvb/UYr7MQoRJH6D9AwmmST1b7EMkhmzZEHSySFIIkHq703trx9pIN1qkgCSrLgf\no/AEkpmkD8dNadyC8qsU2CmCJNVKJ0kHqR0sPwxznZAg3Zom6XazvufFuCJVkGo5+3aUVgep\nj+fRIDmapD6dD50jKQNmQflVCu0UQFKrpZL0oq4cOBqeliBBqknaC0iVhp87B4kmSQ0mK0lp\ngNR1Vsk6K42ScI+8TJHwIOwwyq2BdKj7ftsHSfpl4eh6jPeCLAmSmyP1zGefIBH1kknqJi0T\nKerOlQoBpCFyhQuzAEnRNeY7G0SO5oFkIokGSciDV1B+lQI77SDJJL2ICMkXZMUJIsXIHbp4\nfZNU7gkkS4N0zL6M66ICaUTPThuLM5NkAGnIhFdQfpUCO9u6CxUmdoWiZoUIkjhBpBy5wzxd\n1ZZ7A8l2gmR7s3LSIKnhMw4knaSEQNJniTGT1C8XQBLGFrQnyvsVedsklfsByTrQ8J5Rb0Zq\ntB2QjCRtDKRCBUncgzpJYoV7kKRBOj1yRZAON8sr/MqtgXQ8Wm9teD8Tt6s2ihIkBkfE/Qqm\nzp0RpD4fVkH5VQrrVEES9omr7t1onTTaTUVu/8DF3kCyK+bHKKQgmAeSiSQ3SAWjoPwqBXYW\nzZWxrrLSTqHrLs7IoN3sTUeuANLB8gq/EiD169i58bUgSJzO3WZBonaVCyTt+qsBj9bxCpCY\nigmkMadIJEg0SWaQ1L5dCiAVEkj6fqFBEscWJJCMeAwgHczvwiwBUqdYQJIwGAESo3OngEQ9\ntsMoKL9KYZ0tSHVN6T2jgySNLUggWfBonK8AaVAKXTsRJA5HNEgkSS/Keup9me6C8qsU1lnX\nvd1Xhl0j1516QV+3wIpHD9IBIDXaEUhU506ab8hAkrug/CqFdTJAEnZKvZeMIDnOfG7dXeCF\nve0KD1K/IIau3f1CPGk+DaRurF39NJQrKEgESfQzOeQ1yo2CJO8lE0iuyO1BOgAkSb7mbOie\nflI/TeUKCxLRuaOeyREIkm+uSQukurQvtlMkNkilM3L7u8CtwxK7A+nLV9fOI0giDc8D9jIf\nJMMzORJfw/c+EVdB+VUK6xxA4t//QYPEidz+uSSAVKk/RfI5i1A4kFgciZNsOUh6oR/KJpqk\nPYHEitzaWfftzFecdggSdQ+4Z5DqZw3/+8fWR6uXpyqQKn2M1YuqQ5eU8PtAO7o0+EVeV02F\n6+Kz90/Dj/CloeKZWhOgw+c/6fP55dY2Sf0mRucfewHfqWUSSOO6do/ryeMEkcfSa4uUz2mR\n9DtjpAboRT9t0u+wcRWUX6WwzqZF4jZIxhaJ2wTcuiYJLZKg7EQsG5lGp7hB0h7JoR0qSfGD\n1FR3Lkj8yP3rmqSuNwiQPA42lCI7FpAYNwz4AclJkr5eu1nNUVBFSYM0InKbW4X6JmnnIPkf\nbDgO/3sDqQqOFzZHyozedpIojtS71RwFVbQuSNzBbxKkUZHbgJQDpFodRxfi+b5pF2SFj1hB\nGmS6aVVvkhICacxFAhGkkZH719xz1zVJOwfJpknXkbrHCB13NrhBEkGYBZKNJAtISpOUBkjs\n2xo0kMZHbgtSDpBcCnmvHR+kroFYFCS9bxc9SE1V9Z6d8CfJCFI5JXL//pqbV5smCSCZtRWQ\nLCRtESRi0mUrSFMjtwOpbZJ2DlImSl3Hzo2vWSAxOfIEkkxSoiBJaymQpkdu/VH1CgFSuSJI\nzntBBQrq6JgBkhkl8xOyCYNEcFQSKH0UxazIFUE67B6k8q2aReh+9vYYhUNskEQG5oNkImk7\nIPUciSBRnl51mnNBqie6q0FSL0Hx09wESG8NLYu/+nJpkAwkySCJNrlvlxBI1DSxik24MuYL\npMPeQcraCSKXfkKWDVLfOswEiSZJBEmxKU1SgiBx0pwNUj1jZAwgCUD52J+d+C1S07XLzvq6\n2EDicmR4V7kdJM2WLEgsjvxFbgvSQb1Nb1Im6YJ0b9/qcvR497dNIUDSTqRpkCiSXlSMBttB\n/useOUhtyQWQ7H6PkVs0d9ztHKTy8X7KstOVmgF8TZDEsG7aBiNIkt0GknnaXtImnyQlAZIw\n1uDwAySmor8g63jMhw+S4reBZJi21zAlqdy3SwWkF2m3GOUzcvsnZQESrSRA0jawgURP22tw\nyX07gGQD6QaQHteqa/dOrYoBpKFpeCko6VvYQSKm7TW6pCYpbpCKNUH6A0jCYIOvxygcCgeS\nQJINpEIFyexKDaS6Df1gcRQCpIP8lPqkTNIF6ZJdq2tI1+yir1sRJCKgSZDojawgKdP2WkxS\n3+6DXcVVQcrXAKkladcgVRdiu3/aOnZufE0AqX3IjgLJsJUDJIEks1NtkqIGqQBIAIkUDyTT\nZi6QBtmHJSSQnLHZKQaQnAn6BelP6Ns5nFsF6Zx9rQJ8Mz8AABVhSURBVNW1s02FQICkx7xx\nO48gHdIBqSnuWiDddg7Sz5OiarDhtPxgAwek3AKSeTtfIIlNUhog5cWaIB12DFL5e36C9PZJ\nrVoPJCKaeSCV/kAqVJC4JC0PUqGC5E7QM0h/zSwo4mzHkzLxVRwXUZu6IBsEpNI3SAeABJBm\nKg6Qurnn1Jg3ZOAVJKFJ+mDFp7lKIZ0RgPR3GyZB2SVIa72xLxhIpTeQilRAGjg6sMsZAKQC\nIMUGEhHLuRbzxix2DFK+Pkg3p3ObIDW6X8+/+tLQIBlnuRJD2QSSOQvPILUkMbtMpiqFdEYB\n0t8wCcqOQXpCcySWjUyDowVAKn2BJDVJCYHESTAESDeAVJaXxR81DwdSGSlIWun9gDRwdChW\nBKnsZxPaNUiUogCpG2tQYt6aiV+QGpKYFzqdVRJr4Buk3A1SwMjtZxMCSIrWAkmN5FyIZBZI\n/zyBJDZJM0AyJh8GJGtCQUFqZhMqAJKidEEyx+5MkDgkSVVy5M7aZ4bdJOewPkhlN5tQJCDZ\ngNoZSLkRJHsm/xwBPBKkg9CzZFSxqxIvd95uc4LUnSKtC1IzCUoBkGQFB8kwgakYx1NBYpLE\neHJpKkjc3Fm7jQVSc73LnlDQyJUmuOOnCZAmaDGQeCSNBokR96NzZ+y2VEC6ASRCK4FEhHFe\nyDHvyOSfko4zlI0a+nbjQOLkXYyFk1RXzphAujmdAGmuFgSJ3bmyqS/CmKBnZKzk7k7TpK6Y\nHUiOEAkbueIEd/w0AdIELQkS83TfKg0kN0mTrmIxa2TaW/1Yw9ogNfNyjXsLIECaoDEg5SaQ\nXJl4B+kg3nzuyHziDUq8Ghl2ljDp96ogUXM3rA8SARRA4oTykOioUKbVlUGYe8WWNStNOndW\njeidFRFIzX1CN4AkKDxI1Nz0egwPzzEwQq7SkOj4UFYdBEi2Jzi0NPm526rlmEdzmD3fFSKL\ngHQASJJiBcmZiZAoP5RNm7R9O2kqSTpbGg9e7o6K2UESGqSVQSIeOQdImwDJQRI9FidZ8u6l\n6o6QN+Lhzt1dM/tkS9GBJD4pC5C2AZLrbjc6RcHSg2QlyYaHK3dbulSNtFyF9yIZnMtFrvbI\nOUBaByQpgkmQ3JnIibpDmUhC8DR9O3W6fXMW8x7i4NRIzXZ4dazJuTRIwgN+AGkrINlIMt/t\nJpXDDhIDD0vuqjg1kjMWe3arg1S95hwgyVofJIGjGSBZSDLff62Ww/QmJSYeY5yMGklFFHt2\n8YBURAdSv2CDIGmXPYR48geSGSXL86RiQYY5UKbjMcLJqVFJ9exiAKls+3YAqVecIDEysT8t\nKIeyMRFjQabiMcI5qka5eIoUBUjSk7IAaQ2QxGjyCxKNku1eabEkXkDKVdmQc9eoT1bs2UUA\nUjt5A0DqtS5IuW+QKJTccyE0RTmMIInEQ4Ookvl9gWpNrXdSiT27OEASn5QFSEFBol/NqoSe\nFL0eQFJJYs0XRxTFLg0kEqIGJEeCthoNqccGUik9KQuQ4gSJkwmvnXE5pw43dCAZ+RFBcqZo\nqtFQNqlnFwlIwgN+AGl5kKQICgRSlw3DaS4LAyQ3Rg1IPJQcIOVxgSQ94AeQFgFJOttXOAoF\nEtupgMRskj6YGHUgudO0zZWXW0FaKXLFB/wA0jIgmW4RMIDEysQ7SM/SvLCbpA8mRj1IjDQt\nPeChZxcPSOIDfgApLEgDHv06KXRUjtYASfy7/8Jukl6M4JhAcqcq3IGhrFEapGhA6h/wiw8k\nASi3EgKphyNikIpn0DNByq0gHRqpILnStVybihIk+ZlzgORfFEgtHSpHEYAkNEkvvGtJCh46\nQoNU5zSQlJ5dPCAVAKnTYiA1IcABiZdJGJCGF2NYOaJA0iAaYOKhZAJJbJAiBKl9VBYg+RcN\nUhUDekRGANJQKgEkY8DrHTYNIm0hjyQOSFSN1ovc/lWYACkwSOqNPw6O4gDJSlKugWRufISV\nHJQMIOVazy4akPpXYd4A0pIgETEZHUhDk0TFe8/CCwOi1slGyQZSHitIBUBqFASkIXbqo9+F\nWhcaQujV/w/rm/G97neX3szf/7j+4kNmgygv0ZHTUZJ+v+jNkil98rfsp8r/1/3uYkf6/Sf4\n2wXD9s0CIT3p95/6u10g/B4K2+RXJyHmT5RP/u0un7G8WvlKrXy8eAFIrN/jQBJPbVggmUHx\nAZKYXdcgKeVfF6TbsCcAkmfxunZdTKldO2YmHrt2wmx1uXiWZFYfOy6j0gkUV2g7xPRoxkHt\n2cXTtSv7l8r299zF1rXrv1iUAEjq/Kk0SGIscTMJBlLu5INPUS4MSxCbuEHqy6ODxI/HkJHb\nPJf0CpDWBUlpkFYDqeyzz51N0hiKRJCIRklmyQBSvVERK0j9c0k39nteghbH6LQJIHlziiDZ\nSWpwmHKvHdkoyZU3cSQ3SFGB1D9OwX49RdjimJw2JQ1SF1lxgDRMjW8nqWtVjCBpdZSchsbM\nAVIeOUgVSXmCIPX7cysgDQsdr+sWFRIkuukQOSBBoqqp3iduQUnfTbncIMUJUkvS6+3GnMM4\ncHEMTuqQd99SAEl9D1+sIJV9KBtJEk+ONJBIiJo0zanICZDzQFANUpQgVZdleVOvBi6Oyykc\n7+57yiAJHMUIkk6SDAD/cT39WVoTSvJuGryxgySSxJkxMnRxHM7hcMcP0r9GL5U+aNXrqrN2\nxfJvLSlFe+lu7TmIP15kGepmqm4vOjUhRdE4LJUL3ESG8MW4gO/8N8n5BOij7dw9S7l6cRxO\n4Wh3v6IFqf1ktkjin29+Jr5bpOGFlkrLMUhoPcbPEk62SlKzRD+aMbR7So2ayChjaAJu/W3g\nBWNaoeDFcTqbgibQIrWfLpAOcYJkIEmKcMe0jwRIKkt6ujpI8g5SasSPx/CR277n5bUGyTWJ\nQ/jiOJy31EDq/pzS47rU7FdrglRqIIkRr4T4xLm/rSiRIFENUnwg/bX33L3WIDkemV2gODZn\nPSiSFEiH8SDx6QgJknuGoOmT6FOQmkYCTQ1ShCC1twpVt2Dcbo4n/ZYojtHZli81kA42kPSe\n3boglVyQdDz4IGmJDyypIB3kPzRqjfjxuETktiA1V2br4btRxWlPr9oYly/ueiw4mUn0IJVu\nkJQGKRaQrCjRePBBMqL0oi+lOYoRpKrH9Ho4tIN3N9vt4GISXUhrIAlD6d4Krmbymj5IXajE\nBVL5oRbRgNHs9yMZWyV5gbh7+kKOx2MRkOpT+KbYVLNCFEflRgOpHbjwVHA1k9fXxECiTiiI\nnksdLBGBRJJkx4MPkp66OPIgcaVyFC1IdaTmHUkyHmpxjNxQC7wUXMskTwukgxUkOTgjA0kP\ndjcefJAsKEnNU79rekULUoNHU/hXCx4ObrQF/LeZGVdoab7mT5AOqYCUm0DqAkeJztVB+meJ\ndSYefJA0lF6oCbySA6l4bavwasKD4ubZPORVbL+SIHFfeUEVh8q1zu1w2AZIRM9ufZAs70Zn\n4zHGqYAk7hcTRxGD1I2K9X8O3CDVIf16kNUDVQhOfnFcrdxrpSanxEESGySpZ7c7kCSUiFuE\nhF3TK2aQqvm5bt3wXdu0GkHKhYgm9JpLIN0Yc35xBjByIdO/hECqaEkIJC5J3kDS58qjQBIL\nGTVIZX8OIiLx2jQ8yqcRIbFpUmiwFYdz3lX3H/tiJQNSToIk9WDEFVsASUuTg5wJJDLVuEHq\nb8HJ1camA4hqfF77cyTV8fwt00DnWtDcyAuaTITipA1SbmqQogCJSZKGhy1NToIUSHTqsYPU\nXbC52XtuPSm5dErUt1miIVfxMIwpWEGSilNnksydDSNBKpMFyZ07J03zC2CkBKMH6U9pAkww\nGQbphk15ThZIMpx1uW7p3CLUdeCI65yHaEHikSRUiZk7L02CovRA0sfNRCyUMTkLHiKD1UYj\nQZK6i3JSzwUp3WvXNTw0SNTf3RhAYpE0vGBjRO7MNDXJCSYA0l/9XIUp6B00yAuE5uw154Ek\njWioQ4PtgrRuWu2AEZ+Cy209u7RAmpD7fkAq9ZvbJoF06247aliyOeUGiBwaVK76pgXSQXwu\ne+AoXpA4JE3OfQJISoLJgGS6t4Dipt+UbGfE0yw6zbxvgAwnZMPQoJZrEiDlHUgtM0qDlChI\ns3IfC1K33cgYjwEkGg8REw6DRTPkJnbPtIbHeXG3dVL4pgZSRU2ugKQGTBwgOUianfv+QJpU\nnEI789HbF+fFXRE5QzuYHki5wBHZs0sBJC+5jwCp32Z3IJXtxd0++HPtIq+5Byc3VYZxisRA\n0mecN/XsYgHJiJK/3LkgDVvsESTtNMveg8tNN5CTC7pM0gApt4GkRUw0IJGR7jd3FkiCf58g\n/Q1dPNtF3m6wgT80OGQSP0iFASSqZxcdSGqkh8gdIDEzoQfplB7cmCF1MZOkQaI5igukRZx2\nkETnnkH6cz9l5OzJCQvkTBIAyfD+O4AkygKS5Ns5SH/WuyVGgKRlkghI+itZKY52DFItEiTZ\nApBK890SXJCoTFIFiTxD2jtIlRx3lAOkoYs3DSRDJimAVOggDV09Imx2DZLwEAe1eq3I9eP0\nWxzDmc+k4pQJgXRwcQSQ3M5VIzcykHwWp0wFpFwFKQdII5zbi9zIilMmAVLXJB1YHAEkXduL\n3MiKUyYDkjCT9QEgjXVuL3IjK06ZDkj961hNHAEks3N7kRtZcco0QGqnx1E5ohskgKRre5Eb\nWXHKpEBSTpUAEtu5vciNrDhlaiBVKDk4Aki6the5kRWnTASkwj6pNUByOLcXuZEVp0wZJHqo\nASAJ2m7kRlaccjMgCVsCpF7bjdzIilOmAlJhm9QaIJmc243cyIpTbgUkcUuA1Gu7kRtZccp0\nQTI2SABp0HYjN7LilMmApLxmwcLR5kAiDpojzR1EbmTFKVMFqYgDJGE/ektTS5s4aI40dxC5\nkRWnTBSkIn2QPBxFYzlDhMqIoArmjLc4ZTogfQgvSrVytABItv25IEjhgyreyI2sOGVSIHUo\nqRylBFJaQRVv5EZWnDIxkGqUHBwBpMQy2UJxyuRAoqRsuSZIwo6llVZQxRu5kRWn3AJI6pbh\nQOLvYE0hj+KI4x1jJlsoTgmQxjj5O1hTyKM44njHmMkWilNuACRtyxhAWvQohitOvJEbWXHK\n9EHStwRIaWWyheKUyYNEbOkfpMiPYrji7LHOk4pTxg9SCZBcC8IVZ491nlScMgGQShtI1Jbe\nQErlKIYrzh7rPKk4ZQoglWaQyC0BknnB/yqFziSyOi9RnDJtkOgtZ4OU2lF0O//XASR9/i/h\nyI2sOGUSINFv0AJIfCdAClycMg2QShokw5ZTQUr2KLqdAClwcdSAnajgIJXOl6MKGgtS8kfR\n7QRIgYujBew0rQKScUuA1C/oeDGC1K/YTp1XKY4WsNMUHqRSBcmyJRekzRxFsxMgLVQcPWAn\naQGQ5LcMW7d0gLS9o2h28kFy9vWSqfMqxSECdooWBsmxpQGk7R5FsxMgLVQcImCnaAmQZk/O\ns92jaHaOB8m4IJk6r1IcImCnKGaQdnAUtRUOPIwrbAv4Z1MAabJiBGlHR1FbERQkZhO1r0NA\nBOwUAaTxmdDx6Kc4i4DkaKISOAQei0ME7BTFBFLkR9ERj/OKw8SD7wRI3AURgXR8Svg5EqRU\njqI7xo1pGiOXT0MQkLovqRyCEMWRA1YJZb58gHTs/2tUl0soZ6CxOL5z3lFcN8aXzSTSQxCy\nOKUIkhrKfAUHSShvo3V2Gz+o+k0ji/FFMwkZuXwnQCrZ/fFAu22FyOU7EwDJ6PRwCEbEeDBn\n9CD9B0FpKi6QZsrH8MkiiSLNbaYJkBZOFGluM02AtHCiSHObaQKkhRNFmttMcysgQdCqWhek\n6ZeDISgurXpnAwTtXgAJgjwIIEGQB60LUtcj7XumR2KZj0S7M8ipiYYo6FJpxlz3oXBh0pxb\nd75WBamr6VH4oS/zkGi1J6WfERR0qTRjrntfuEBpzq37CMUF0pE6bh4SFY6Tn2DyUNCl0oy5\n7hRAPtOcW/cRigKk/odPkMREy3LeDjWk6esvaIAALdOoO1Fnn2mWM+s+QgBpRpoAaW6aAMmL\nxPodhX/Nl6nnnESi/cfERA1ptsvnptlsL1a+qb+fNEtpn3pLczhgM9IU6z33wFNp9h87GWzo\nPtRYmt0iUcE0968yFUw+/tKrweQrzX6J3zQHtqanKTYivlokKc1STG4fLZJcf08gqfHvI+hD\npKmD5DFNKaW40jw28guSmGY58xiNUCQgDUHqESQh8r0FvZimTtaUNMVPT8Ekf3qqu5ymn7p3\nHz7rLqZZzqz7CMUBknpYvIAkceQJJJ0jXwHqpfJUmqWnuktpeqq7nLT6Y36a5cy6j1AUdza0\nrfFRXDb7zgYp0VJJ3EOaXRfC0x0DfipPpemp7mKa3upeio2c1zsb/NSdL9xrB0EeBJAgyIMA\nEgR5EECCIA8CSBDkQQAJgjwIIEGQBwEkCPIggARBHgSQIMiDABIEeRBAgiAPAkgQ5EEACYI8\nCCBBkAcBJAjyIIAEQR4EkCDIgwASBHkQQIIgDwJIEORBAAmCPAggQZAHASQI8iCABEEeBJAg\nyIMAEgR5EECCIA8CSBDkQQBpQWW1jte7vPjT+qKE+znLTu3m9HIitQzHdWFhhy+orNO3sti2\n0bHagPANy9VMnElC/oUdvqCa8L5fsuNDX2zfaP5yKKiw1xdUF+OX7P35/89b1c1r2qnnz8cl\nyy4DYPfq571txYbNs+z+1m+ViVsJy9ucuiSGjaBgAkgLqiPiNzuX5XfTy7t2INVdtf6k51H/\nfLZcGkjHYatM2KrZ4E0AqU9i2AgKJoC0oPpeV/XllH1VSHXNx3sV59fss3VcK9bO2VXsqjXe\n86P8zI7aVtfsUv70qVX/iUl0G0GhBJAWlATSs+v1/X7uQ/9UL3o2KY1O2bNPdq/aGhWke/9N\n2uqUPYS0G1KHJO5i7lAAYecuKBmkc5YJJzTdgJ5k7deW8i9iK3lLwgqQggo7d0F1ofxTNSGX\n7PT5fQdIGxF27oLqQvmtOqlphuqUTlovc9dO+jZsZe/aiblDIYSdu6CG60j1j5/yMZwjXatR\nga9qfKCWebBB+jZsVX37FVoqPQmAFFLYuQuqv7Php6zifOiVHbvB6uy3tYpj18PmGkjDVvdu\nIDzrhvTUJABSSGHnLqiGnNO16YRdsuxcD1g3A9P3+nfvFa6mDpvrJz7DVr/nZoNhbFxJAiCF\nFHYuBHkQQIIgDwJIEORBAAmCPAggQZAHASQI8iCABEEeBJAgyIMAEgR5EECCIA8CSBDkQQAJ\ngjzo/20ppMQ2xnWTAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 300,
       "width": 420
      },
      "text/plain": {
       "height": 300,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "scaling_parameter=max(figure_data$domestic_backproj+figure_data$imported_backproj)/5.0\n",
    "options(repr.plot.width=7,repr.plot.height=5)\n",
    "\n",
    "figure_data %>% \n",
    "    select(date,imported_backproj,domestic_backproj) %>% \n",
    "    gather(exp_type,value,-date) %>%\n",
    "    ggplot(aes(x=date)) + \n",
    "        geom_bar(aes(y=value, fill=exp_type, group=exp_type), stat='identity', width=0.7) +\n",
    "        geom_line(data = figure_data, aes(y=Rt_median*scaling_parameter),color=\"#1380A1\",size=1) +\n",
    "        geom_ribbon(data = figure_data, \n",
    "                    aes(ymax=Rt_upper*scaling_parameter, ymin=Rt_lower*scaling_parameter), \n",
    "                    fill=\"#1380A1\", alpha = 0.25) +\n",
    "        geom_ribbon(data = figure_data, \n",
    "                    aes(ymax=Rt_IQR_upper*scaling_parameter, ymin=Rt_IQR_lower*scaling_parameter), \n",
    "                    fill=\"#1380A1\", alpha = 0.5) +\n",
    "        scale_fill_manual(values=c(\"#FAAB18\",\"gray35\")) +\n",
    "        ggtitle(\"Entire Japan excluding unknown cases\") +\n",
    "        labs(x=\"\\nDate of infection\", y=\"Incidence\\n\") +\n",
    "        theme(text = element_text(size=12, family=\"sans\",color=\"black\"),\n",
    "              axis.text = element_text(size=10, family=\"sans\",color=\"black\"),\n",
    "              panel.grid.major = element_blank(), panel.grid.minor = element_blank(),\n",
    "             legend.position = \"none\") +\n",
    "        guides(fill=F) +\n",
    "        theme_bw() +\n",
    "        scale_x_date(date_labels=\"%m/%d\",date_breaks  =\"10 day\", limits = range,  expand = c(0, 0)) +\n",
    "        scale_y_continuous(limit=c(0,620), expand = c(0, 0),\n",
    "                           sec.axis = sec_axis(~./(scaling_parameter), breaks=c(0,2,4,6,8,10), name = \"Effective reproduction number\\n\")) +\n",
    "        geom_hline(yintercept=1*scaling_parameter, linetype=\"dashed\", color = \"#1380A1\", size =0.7)"
   ]
  },
  {
   "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": "4.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}