{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“replacing previous import ‘vctrs::data_frame’ by ‘tibble::data_frame’ when loading ‘dplyr’”\n",
      "This is surveillance 1.18.0. For overview type ‘help(surveillance)’.\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "'R version 4.0.2 (2020-06-22)'"
      ],
      "text/latex": [
       "'R version 4.0.2 (2020-06-22)'"
      ],
      "text/markdown": [
       "'R version 4.0.2 (2020-06-22)'"
      ],
      "text/plain": [
       "[1] \"R version 4.0.2 (2020-06-22)\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "libraries = c(\"dplyr\",\"magrittr\",\"tidyr\",\"ggplot2\",\"readxl\",\"RColorBrewer\",\"surveillance\")\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",
    "theme_set(theme_bw())\n",
    "version$version.string"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Revision note\n",
    "\n",
    "$\\circ$ typo in the function of generation interval (shape parameter) was corrected. (15 May 2020)\n",
    "\n",
    "$\\circ$ assuming the proportion of asymptomatic cases among all cases is constant, the asymptomatic cases were included in the estimation of $R_t$ including all cases with unknown illness onset date. (16 May 2020)\n",
    "\n",
    "$\\circ$ the convolution between the incubation period distribution and time distribution from the illness onset to reporting was corrected. (22 Nov 2021)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "For the analysis, all confirmed COVID-19 cases who diagnosed in Japan were categorized into two types of case, imported case and domestic transmission case. The imported cases were defined as (i) those who have international travel history, (ii) who embarked on the “Diamond Princess” cruise ship as a passenger or member of disaster assistance team (i.e., DMAT or PMAT) within two weeks before the illness onset. \n",
    "\n",
    "For the estimation of $R_t$, the symptomatic confirmed cases in Japan with either or both the date of illness onset or the date of laboratory confirmation were used. Since the cases from the chartered flights have not involved in the domestic transmission dynamics in Japan, those cases were excluded in the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## The code was written by Sung-mok Jung\n",
    "\n",
    "## load the data\n",
    "datestar = as.Date(\"2020-05-10\")\n",
    "read.csv(\"../data/JapaneseDataCOVID19 (\"%&%format(datestar,\"%y%m%d\")%&%\").csv\", encoding=\"UTF-8\") -> master_df \n",
    "colnames(master_df) <- c(\"type\",\"onset\", \"labconf\", \"reported\", \"asym\")\n",
    "master_df %<>% mutate(labconf = as.Date(labconf), onset = as.Date(onset), reported = as.Date(reported))\n",
    "\n",
    "\n",
    "## data for cases with the illness onset\n",
    "master_df %>% filter(!is.na(onset)) -> onset_Japan\n",
    "onset_Japan %>% group_by(onset) %>% count(type) %>% spread(type, n) %>% as.data.frame() -> master_dff\n",
    "master_dff[is.na(master_dff)] <- 0 \n",
    "master_dff$onset <- as.Date(master_dff$onset)\n",
    "    \n",
    "as.numeric(datestar-1-min(master_dff$onset)) -> time.diff\n",
    "ttime <- as.data.frame(c(0:time.diff))\n",
    "colnames(ttime) <- c('t')\n",
    "ttime %<>% mutate(onset = as.Date(min(onset_Japan$onset))+t)\n",
    "    \n",
    "merge(ttime, master_dff, by='onset', all.x=TRUE) -> master_df_final\n",
    "master_df_final[is.na(master_df_final)] <- 0\n",
    "    \n",
    "master_df_final %<>% mutate(total = domestic+imported) %>% dplyr::select(onset, domestic, imported, total) %>%\n",
    "rename(domestic_org = domestic, imported_org = imported)\n",
    "master_df_final -> entire_Japan\n",
    "\n",
    "\n",
    "## data for cases with the unknown illness onset (excluding the asymptomatic cases)\n",
    "master_df %>% filter(is.na(onset)) %>% filter(!is.na(labconf)) -> unknown_df\n",
    "unknown_df %>% group_by(labconf) %>% count(type) %>% spread(type, n) %>% as.data.frame() -> unknown_dff\n",
    "unknown_dff[is.na(unknown_dff)] <- 0 \n",
    "unknown_dff$labconf <- as.Date(unknown_dff$labconf)\n",
    "    \n",
    "as.numeric(datestar-1-min(unknown_dff$labconf)) -> time.diff\n",
    "ttime <- as.data.frame(c(0:time.diff))\n",
    "colnames(ttime) <- c('t')\n",
    "ttime %<>% mutate(labconf = as.Date(min(unknown_df$labconf))+t)\n",
    "merge(ttime, unknown_dff, by='labconf', all.x=TRUE) -> unknown_df_final\n",
    "unknown_df_final[is.na(unknown_df_final)] <- 0\n",
    "    \n",
    "unknown_df_final %<>% mutate(total=domestic+imported) %>% dplyr::select(labconf, domestic, imported, total)\n",
    "\n",
    "temp_start <- matrix(NA, ncol=1, nrow=as.numeric(min(unknown_df_final$labconf)-as.Date(\"2020-01-03\")))\n",
    "temp_start %<>% as.data.frame() %>% mutate(labconf = as.Date(as.Date(\"2020-01-03\"):(min(unknown_df_final$labconf)-1), origin=\"1970-01-01\"),\n",
    "                                           domestic = 0, imported = 0, total= 0) %>% dplyr::select(-V1)\n",
    "rbind(temp_start, unknown_df_final) -> unknown_Japan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributions\n",
    "\n",
    "### 1. Generation time\n",
    "In the present analysis, serial interval was assumed to be identical with generation time and the value of doubly interval censored serial interval was also adopted from the previous study (Nishiura et, al., 2020, IJID). The Weibull distribution was fitted to 28 known infector-infectee pairs (mean = 4.8  days, standard deviation = 2.3 days).\n",
    "\n",
    "\n",
    "### 2. Incubation period\n",
    "We accounted for right truncation and modeled and used a lognormal distribution with parameters adopted from an earlier study (Linton et, al., 2020, JCM), and the estimated mean and SD were 5.6 days and 3.9 days, respectively.\n",
    "\n",
    "\n",
    "### 3. Right-truncated reporting delay\n",
    "We fitted the right-truncated time delay distribution from the illness onset to reporting to the Weibull distribution. We define the likelihood as follows: \n",
    "\n",
    "$$\n",
    "L(\\theta|S_k,O_k,T)=\\prod_{k=1}^{K}\\frac{h(S_k-O_k  |  \\theta)}{H(T-O_k  |  \\theta)}\n",
    "$$\n",
    "\n",
    "where $S_k$ and $O_k$ are the reporting date and the date of illness onset of case $k$, respectively and $T$ is the latest calendar day of observation. $h(∙)$ and $H(∙)$ are the probability density function (PDF) and cumulative density function (CDF) of Weibull distribution and $\\theta$ is the parameters of the distribution (i.e., shape and scale). The likelihood is maximized to determine the best-fit parameters used in the estimation of the reproduction number ($R_t$).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Serial interval (Nishiura, et al, 2020)\n",
    "gi_fit = list(shape=2.305, scale=5.452)\n",
    "generation <- function(t){pweibull(t, shape = gi_fit$shape, scale = gi_fit$scale) - \n",
    "                          pweibull(t-1, shape = gi_fit$shape, scale = gi_fit$scale)}\n",
    "\n",
    "\n",
    "## incubation period (Linton et al, 2020)\n",
    "inc_fit = list(meanlog=1.519, sdlog=0.615)\n",
    "incubation <- function(t){plnorm(t, inc_fit$meanlog, inc_fit$sdlog) - plnorm(t-1, inc_fit$meanlog, inc_fit$sdlog)}\n",
    "\n",
    "\n",
    "## right-truncated reporting delay from illness onset to laboratory confirmation (estimated values using the MHLW data)\n",
    "rep_fit = list(shape=1.741, scale=8.573)\n",
    "onsettolabconf <- function(t){pweibull(t,  shape=rep_fit$shape, scale=rep_fit$scale) - \n",
    "                              pweibull(t-1, shape=rep_fit$shape, scale=rep_fit$scale)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Back-projection from the laboratory confirmation date to illness onset\n",
    "\n",
    "To estimate the effective reproduction number ($R_t$) by time of infection, we back-projected the cases by time of infection.\n",
    "\n",
    "First, since there are some reported cases with unknown illness onset date, we back-projected those cases from laboratory confirmation date to illness onset date, using the PDF of time delay from the illness onset to laboratory confirmation. $R$ package ‘$surveillance$’ was used for the back-projection and within the package, the non-parametric back-projection algorithm based on Expectation-Maximization-Smoothing (EMS) algorithm was conducted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## adding extra 10 days for the stability of back-projection procedure\n",
    "temp_lastdays <- matrix(NA, ncol=1, nrow=10)\n",
    "temp_lastdays %<>% as.data.frame() %>% mutate(labconf=as.Date((max(unknown_Japan$labconf)+1):(max(unknown_Japan$labconf)+10), origin=\"1970-01-01\"),\n",
    "                                              domestic=0, imported=0, total=0) %>% dplyr::select(-V1)\n",
    "rbind(unknown_Japan, temp_lastdays) -> unknown_Japan\n",
    "unknown_Japan %<>% mutate(time_onset = 1:nrow(unknown_Japan))\n",
    "\n",
    "\n",
    "## time delay from the illness onset to laboratory confirmation\n",
    "K = nrow(unknown_Japan)\n",
    "\n",
    "report_probability = pweibull(1:K, shape=rep_fit$shape, scale=rep_fit$scale) - pweibull(1:K-1, shape=rep_fit$shape, scale=rep_fit$scale)\n",
    "report_pmf = c(0,report_probability[1:21])\n",
    "\n",
    "\n",
    "## back-projecton of imported cases\n",
    "sts = new(\"sts\", epoch=unknown_Japan$time_onset, observed=unknown_Japan$imported)\n",
    "bpnp.control = list(k = 2, eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                    Tmark = nrow(sts), B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, eq3a.method = c(\"R\",\"C\"))\n",
    "sts_bp = backprojNP(sts, incu.pmf=report_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "unknown_Japan$imported_backproj = upperbound(sts_bp)\n",
    "\n",
    "\n",
    "## back-projecton of domestic cases\n",
    "sts = new(\"sts\", epoch=unknown_Japan$time_onset, observed=unknown_Japan$domestic)\n",
    "bpnp.control = list(k = 2, eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                    Tmark = nrow(sts), B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, \n",
    "                    eq3a.method = c(\"R\",\"C\"))\n",
    "sts_bp = backprojNP(sts, incu.pmf=report_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "unknown_Japan$domestic_backproj = upperbound(sts_bp)\n",
    "\n",
    "\n",
    "## normalizing the back-projected cases\n",
    "unknown_Japan$imported_backproj[unknown_Japan$imported_backproj<=0.01] <- 0\n",
    "unknown_Japan$domestic_backproj[unknown_Japan$domestic_backproj<=0.01] <- 0\n",
    "\n",
    "unknown_Japan %>% mutate(imported_normal = imported_backproj/sum(imported_backproj)*sum(imported),\n",
    "                         domestic_normal = domestic_backproj/sum(domestic_backproj)*sum(domestic),\n",
    "                         total = imported_normal+domestic_normal,\n",
    "                         time_onset=0:(nrow(unknown_Japan)-1)) %>%\n",
    "dplyr::select(time_onset, labconf, imported_normal, domestic_normal, total) %>% \n",
    "rename(t = time_onset, onset = labconf, imported_backproj=imported_normal, domestic_backproj = domestic_normal) -> dt.backproj_onset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## merge the back-projected data with cases whose date of illness onset was available\n",
    "merge(entire_Japan, dt.backproj_onset, by=c('onset'), all.y=TRUE) %>% \n",
    "mutate(domestic = domestic_org+domestic_backproj, imported = imported_org + imported_backproj,\n",
    "       total = domestic + imported) %>% dplyr::select(onset, domestic, imported ,total) -> df_onset\n",
    "df_onset[is.na(df_onset)] <- 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Back-projection from the illness onset to time of infection\n",
    "\n",
    "The epidemic curve by date of illness onset was back-projected once again to the epidemic curve by the time of infection, using the PDF of incubation and the same back-projection method was used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## adding extra 10 days for the stability of back-projection procedure\n",
    "temp_lastdays <- matrix(NA, ncol=1, nrow=10)\n",
    "temp_lastdays %<>% as.data.frame() %>% mutate(onset=as.Date((max(df_onset$onset)+1):(max(df_onset$onset)+10), origin=\"1970-01-01\"),\n",
    "                                              domestic=0, imported=0, total=0) %>% dplyr::select(-V1)\n",
    "rbind(df_onset, temp_lastdays) -> df_onset\n",
    "df_onset %<>% mutate(time_onset = 1:nrow(df_onset))\n",
    "\n",
    "\n",
    "## incubation period\n",
    "K = nrow(df_onset)\n",
    "incubation_probability = plnorm(1:K, inc_fit$meanlog, inc_fit$sdlog) - plnorm(1:K-1, inc_fit$meanlog, inc_fit$sdlog)\n",
    "inc_pmf = c(0,incubation_probability[1:21])\n",
    "\n",
    "\n",
    "## back-projecton of imported cases\n",
    "sts = new(\"sts\", epoch=df_onset$time_onset, observed=df_onset$imported)\n",
    "bpnp.control = list(k = 2, eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                    Tmark = nrow(sts), B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, eq3a.method = c(\"R\",\"C\"))\n",
    "sts_bp = backprojNP(sts, incu.pmf=inc_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "df_onset$imported_backproj = upperbound(sts_bp)\n",
    "\n",
    "\n",
    "## back-projecton of domestic cases\n",
    "sts = new(\"sts\", epoch=df_onset$time_onset, observed=df_onset$domestic)\n",
    "bpnp.control = list(k = 2, eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                    Tmark = nrow(sts), 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_onset$domestic_backproj = upperbound(sts_bp)\n",
    "\n",
    "\n",
    "## normalizing the back-projected cases\n",
    "df_onset$imported_backproj[df_onset$imported_backproj<=0.01] <- 0\n",
    "df_onset$domestic_backproj[df_onset$domestic_backproj<=0.01] <- 0\n",
    "\n",
    "df_onset %>% mutate(imported_normal = imported_backproj/sum(imported_backproj)*sum(imported),\n",
    "                    domestic_normal = domestic_backproj/sum(domestic_backproj)*sum(domestic),\n",
    "                    total = imported_normal+domestic_normal,\n",
    "                    time_onset=0:(nrow(df_onset)-1)) %>%\n",
    "dplyr::select(time_onset, onset, imported_normal, domestic_normal, total) %>% \n",
    "rename(t = time_onset, imported=imported_normal, domestic = domestic_normal) %>% filter(onset <= (datestar-1))-> dt.backproj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Renewal equation\n",
    "\n",
    "Let $C_{domestic} (t)$ be the number of newly domestic infected cases (i.e., defined as confirmed COVID-19 cases who were infected within the Japan) on day $t$. Then, the expected number of domestic infected case is modelled as:\n",
    "$$\n",
    "E(C_{domestic} (t))=R_t \\sum_{\\tau=1}^{t-1}C_{total} (t-\\tau)g(\\tau)  \\frac{F(T-t)}{F(T-t+\\tau)},\n",
    "$$\n",
    "where $C_{total} (t)$ is the number of total infected COVID-19 cases (i.e., sum of domestic infections and imported infections) on day $t$, $R_t$ is the reproduction number on day $t$ and $g(∙)$ is the PDF of generation interval distribution. In addition, to account for right-censoring, i.e., cases that have not yet been reported but already infected, $F(∙)$, the CDF of time delay from infection to reporting (i.e., the convolution of PDF of incubation period and PDF of right-truncated reporting delay from the illness onset to reporting) was considered. \n",
    "\n",
    "To estimate parameters governing $C_{total} (t-\\tau)g(\\tau)\\frac{F(T-t)}{F(T-t+\\tau)}$, we assumed that $C_{domestic} (t)$ follows a Poisson distribution and we arrived at the likelihood function of the form:\n",
    "$$\n",
    "L(R_t;C_{domestic} (t))= \\prod_{t=1}^{T}\\frac{\\exp⁡(-E(C_{domestic} (t)) E(C_{domestic}  (t))^{C_{domestic}(t)}}{C_{domestic}  (t)!}.\n",
    "$$\n",
    "\n",
    "The 95% confidence intervals of each estimates were derived using the profile likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "today = as.numeric(datestar-as.Date(min(dt.backproj$onset)))\n",
    "\n",
    "## start estimation from 1 Feb 2020, as there is a huge uncertainty in the early phage of epidemic\n",
    "Start.T = 29\n",
    "\n",
    "\n",
    "## reporting delay distribution from infection to lab confirmation\n",
    "infectiontoreport <- function(t, tau){sum(convolve(onsettolabconf(1:1000),\n",
    "                                                   rev(incubation(1:1000)),type = c(\"open\"))[1:(today-t+tau)])}\n",
    "\n",
    "## pre-calculation for reporting delay to reduce the calculation time\n",
    "precal <- matrix(0, nrow=(max(dt.backproj$t)-Start.T+2+1), ncol=(max(dt.backproj$t)-1))\n",
    "for (m in (Start.T-2):max(dt.backproj$t)){\n",
    "    for (n in 1:(m-1)){precal[m-Start.T+2+1,n] <- infectiontoreport(m,n)}}\n",
    "\n",
    "delay_precalculation <- function(t){sum(convolve(onsettolabconf(1:1000),\n",
    "                                                 rev(incubation(1:1000)),type = c(\"open\"))[1:(today-t)])}\n",
    "dt.backproj %<>% rowwise %>%\n",
    "mutate(imported_delay = imported/delay_precalculation(t),\n",
    "       domestic_delay = domestic/delay_precalculation(t)) %>% mutate(total=imported_delay+domestic_delay)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "## estimating Rt\n",
    "est.t <- list()\n",
    "est.CI <- list()\n",
    "\n",
    "for (TT in Start.T:max(dt.backproj$t)){  \n",
    "    \n",
    "    dt.backproj %>% filter(t <= TT) -> dt.backproj.T\n",
    "    \n",
    "    llk <- function(param){\n",
    " \n",
    "        llk <- 0; Cs <- 0  \n",
    "        t=TT; R = param\n",
    "\n",
    "        Css <- rep(0, t)    \n",
    "\n",
    "        for (tau in 1:t-1){Css[tau] = (dt.backproj.T$total[t-tau+1])*generation(tau)/precal[t-Start.T+3,tau]}\n",
    "\n",
    "        Cs = sum(Css)*R \n",
    "        Cs[Cs<=0] <- 1e-5\n",
    "\n",
    "        return(-(-Cs+dt.backproj.T$domestic_delay[t+1]*log(Cs)-lgamma(dt.backproj.T$domestic_delay[t+1]+1))) \n",
    "    }  \n",
    "\n",
    "    param0 = c(0.7)\n",
    "    opt_est <- optim(param0, fn=llk, method=c(\"L-BFGS-B\"), lower=c(0), control=list(maxit=10000))\n",
    "    opt_est$par -> est.t[[TT]]\n",
    "\n",
    "    ## 95% confidence intervals using the profile likelihood\n",
    "    ci_pro <- matrix(NA, ncol=2, nrow=1)\n",
    "\n",
    "    CI <- function(par_CI){return(2*(-llk(par_CI)+opt_est$value))} #llk\n",
    "\n",
    "    par_CI <- seq(0, 10, by = 0.01)\n",
    "    \n",
    "    logLik <- sapply(par_CI, FUN = CI)\n",
    "    as.data.frame(par_CI) -> par_CI; as.data.frame(logLik) -> logLik\n",
    "\n",
    "\n",
    "    cbind(par_CI, logLik) -> data_CI\n",
    "\n",
    "    data_CI$logLik[data_CI$logLik<(max(data_CI$logLik)-3.84)] <- NA\n",
    "    data_CI %<>% na.omit()\n",
    "\n",
    "    min(data_CI$par_CI) -> ci_pro[1,1]\n",
    "    max(data_CI$par_CI) -> ci_pro[1,2]\n",
    "    \n",
    "    as.data.frame(ci_pro) -> ci_pro\n",
    "    colnames(ci_pro) <- c(\"lower\",\"upper\")    \n",
    "    ci_pro -> est.CI[[TT]]\n",
    "}\n",
    "\n",
    "matrix(unlist(est.t),ncol=1,byrow=T) -> est.t\n",
    "matrix(unlist(est.CI),ncol=2,byrow=T) -> est.CI\n",
    "\n",
    "cbind(est.t, est.CI) -> est"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figures for the estimated Rt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "est %>% as.data.frame() %>% mutate(t = Start.T:max(dt.backproj$t)) -> result\n",
    "merge(dt.backproj, result, by='t', all.x=TRUE) -> result\n",
    "colnames(result) <- c(\"t\",\"onset\",\"imported\",\"domestic\",\"total\",\"imported_delay\",\"domestic_delay\",\"Rt\",\"lower\",\"upper\")\n",
    "\n",
    "result %>% filter(onset<=as.Date(datestar-(17))) -> result_withcut\n",
    "\n",
    "result_withcut %>% dplyr::select(onset, domestic, imported) %>% gather(onset) -> figure_data\n",
    "result_withcut %>% dplyr::select(onset) -> TIME\n",
    "rbind(TIME, TIME) -> TIME\n",
    "cbind(figure_data, TIME) -> figure_data\n",
    "colnames(figure_data) <- c(\"subject\",\"value\",\"onset\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAAJYCAMAAABvmDbGAAAAQlBMVEUAAAATgKEzMzNNTU1Z\nWVloaGh8fHyMjIyampqemk+hzNmnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD6qxj////s\n5NGUAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2di5aquhJF8XjVVvu1Vf7/V2+j\nQN6hklRCImuOcXqrhCKEmodXCF0PAEimW7sCALwDEAkABiASAAxAJAAYgEgAMACRAGAAIgHA\nAEQCgAGIBAADEAkABnKK1MkYU7/3YyFiLNaqAcDLeiKNP0Ek8A7kFSl+alppAAoDkQBgoKRI\nf99vp677+NfPh31joa67H7rT34f7dd/tP27OWN9/878CPH/6PojS2qR5SQAUoKxIt/1Tn3+m\nSH9pf+37cXr344h1HGf7ff30JZXWJoklAVCAsiId/nYf/47dhzR1FOlw7+99v+++7v3t2u3v\n1lhf3fFv93M7dsfXT8PX+0e3t0ySlwRAdkpdtXt9/xr++Td96+d/xr3K16tA/9l9mrH+/hy6\n52HcfQrwlKY/DXPpk+QlAZCdsiLdxt/F3/kcafh8GmtzHxVRYhlfxsO4/vd5dqVNuplzAZCP\nwhcbpH9VkcbPzttO0y/336/TXgkwf7BMgkigEG2J9HtQdnDyNOskiAQKUZlI/li/XXf4/P53\nt4hknwSRQCGqEuk4nvXInKRih/G6+BTgdXH7eY6kT7LXAIBMVCXS53jZQL58sH/KdXtefhiL\nfU0BXhe3T923OcleAwAysaZIv/OX6XLBvjvd+t/rdNFt4KM73of7Q8MF8cNw0/b+N/l5kW/o\nyHAf7iMdLJPsNQAgEyV7f6vpfVJ6Nrym/HRGz4axi8LzmYvv19TPw9PBPxOf0/Y3yyR5SQBk\nZz2Rbh9PPZQpz7523YdypvQsd319/j123el3PPQbOtQdu8OnfZK0JACy03CmwRJQDw0nI0QC\n9dBwMkIkUA8NJyNEAvXQcDJCJFAPSEYAGIBIADAAkQBgACIBwABEAoABiAQAAxAJAAYgEgAM\n5BPpPwAagCndM4qULTIAbEAkABiASAAwAJEAYAAiAcAARAKAAYgEAAMQCQAGIBIADEAkABiA\nSAAwAJEAYAAiAcAARAKAAYgEAAMQCQAGIBIADEAkABiASAAwAJEAYAAiAcAARAKAAYgEAAMQ\nCQAG1hepm96YZ354ApFAA6wu0jDjU5zpP/HhBUQCDbC2SN301/wwApFAA9QhUg+RQNusL9J4\nRgSRQMusLpJyeqSLxDnMPwAZWV2k6S/2SKBlIBIADEAkABiASAAwsLZIlvuwuCEL2mN1kdBF\nCLwD64u0BEQCDQCRAGAAIgHAAEQCK/B4snYtOIFIoCCTPxDJBUQCBCDSEhAJENBFeh+hIBIo\nCERaAiIBAhBpCYgECECkJSASIOASqX2hIBIoCERaAiIBAosiNWsURAIFgUhLQCRAACItAZEA\nAYi0BEQCblRPeojkBiIBN+EiNScURAL5gUhkIBJwA5HIQCTgBiKRgUjADUQiA5GAG4hEBiIB\nNxCJDEQCbiASGYgE3EAkMhAJuIFIZCAScAORyEAk4AYikYFIwE28SM0IBZFAfiASGYgE3EAk\nMhAJuIFIZCAScAORyEAk4AYikYFIwA1EIgORgBuIRAYiARO7J5YJEGkCIgETiBQMRAImECkY\niARMIFIwEAmYQKRgIBIwgUjBQCRgApGCgUjABCIFA5GACUQKBiIBE4gUDEQCJhApGIgETCBS\nMBAJmECkYCASMIFIwUAkYAKRgoFIwAQiBQORgAlECgYiAROIFAxEAiaMItUuFEQC+YBIwUAk\nYAKRgoFIwAQiBQORgAlECgYiAROIFAxEAiYQKRiIBASLOkAkFxAJCCBSNBAJCCBSNBAJCCBS\nNBAJCCBSNBAJCCBSNBAJCCBSNBAJCCBSNBAJCCBSNBAJCCCSh05Fn8pbLQmI1B4QyQNEAlS2\nLpLXhtGd1z8QCXjYuEiGHQoQCVDZtkgd9kiAB4jkASIBKlsU6Y/Xpw7nSICJLYo0f4JIgIst\ni9T1EAkwsWGRuvmPC8UdiAQ8bFkk+31WMhAJCDYs0hOaDVbbIBIQQCQKEAksAJE8KD3tTnd9\nKnvNJiBSe2xdJC9qn9UPfSpvtSQgUntAJAp/h3Y3XLUDHiAShUEiiAQ8QCQKg0SfV/1HltrY\ngEjtAZEo4KodWAAieWjvqt2FsRaAAlkHiNTSVTuIVBqIFEIzV+0uMKkwBUSaf6iMd75qB5FK\nA5FCaeGq3eUCkwoDkUh8n4bLDD/WaRAJQCQax+mKnW1ipSLBpKJAJAJf3f7n79zofuy+LFMh\nEoBIJA7dv+dFhnt3sEytT6TLBSaVBiIReF6os16xe05lq5EORGoHiERgEune7W1T2WqkA5Ha\nASIR+BjOjRo6R7rApOJAJAK3fffqJ3S0Ta1OpAtEKg9EonA7DSK1ch/pApPKA5GSgUgAIjFQ\nsUgwqRgQicTX0Lfh9G2dVptIF4i0AhCJwtRFqImLDReYtAIQicC1O/7+/fN77PSe3wMQCUAk\nEl33ery8jRuyEGkNIBKBuWdQC12ELheYtAIQicCpuw67pPuH9SQJIgGIROK2H3vb7W+WqRAJ\nQCQa9+twcrS/6iNxPYFIACIxUJdIF4i0ChAphAYuNugiwaQyQCQCvwdpfEjjIVmIBCASCdmj\nrtP7gEMkAJFIdN1t+mCbylMdCxEiGR5BpDJAJAL7qUMDRAIOIFIyEAlAJAYqFwkmFQEiEZCu\nNFQ+iL7FI4hUBIhEACKBJSBSELY+QhAJQKQgvo+1X7WDSGsBkaj8fti6NfT1iwSTSgCRSNyu\n+z+LPm1PUUCkLRPoyWKB9xbp++C0qIdImwYiBV61sw36PU1lqIsdiFQ9a4hUm1D0NP05dt3+\nij0SMIBIYWl6/3Qf3VUkktUjiJQRiBR+1e5a/1U7u0gwKR8QKebA6af2+0gQqTQQ6S17NkCk\n0kCkt+z9DZFKA5HC0vT71LXwojGIVBqIFPU2ipNtYopI47xdp394wiUSTMoGRApJ069u/5Pl\nZcyjNt0UZf4QWsMRiFQaiBQ2itC/Z9LfmS9/d0Kf51+hVWgNRyBSaSBS8NsorA/1PadG14Bb\nJJdHECkbEClGJOb3I3U9RGodiBSSph/DuRH/OZJXpP8GAgNCpOJApMDXurwGbmB9h6xyiSHv\nHgkm5QIiBaXp7TSIxHsfqVP/QKQmgUir92yQxiWawkCk5oBIq4sk5oVI7QKRKhKJ7YYsRCoO\nRKpJJK4uQm6RYFImIFIdIvn5Lyz/PR5BpExAJIgEGIBIEAkwAJEgEmAAIkW/jcKcylclDYhU\nPRBpYyLBpDxApIhDu9/jqfDgJxCpdiBSzDkS94N9S0Ck6oFIURcbcGgHVCBSlEhX2+gntYjk\n9Qgi5QEivd/lb79IMCkLa4pUi1ANjGsHkWoHIuUc1+7nNJxMnVzvgSEDkaoHImUc1+44Pq+3\nTzUJIlVLpCeLBd5bpLBx7b66430Q6av7SKofRKoYiCTINa7dvruLOZJgFQkmcQKRBLnGtZuN\nYxApJPshUkEgkiDXuHaHcY/0z3ogGAJEqhaIJMg1rt14jvSz970KnUSQSEseQSROIJIg27h2\np/FiuVW7ECBStUAkQb6eDT/P27ffyZEhUrVAJEELXYQ4RYJJjEAkAU2kTsVWgrdaEhCpWiCS\nIJtI9+twlXx/tT0FGAREqhaIJAjpa7f//fv7s/+0TdRFuu3Hu0gMXYQgUqVAJAFdpNPQRagf\n7gxdLVN1kY7dx7Avul/tfVwDgEjVApEEgT0blA/K1KDSITCLBJP4gEiCwE6r/bBHonQRGvra\nDdyLikTwCCLxAZEEgY9RDOdIpEO7a3cczqh+j9bSIUCkaoFIgoCLDd6+CvbnkUr3bIBIRYFI\ngpAbsj8f7r4K5gHc88H0Y2pPO4hUMRBJ8F49GygiwSQ2IJKgCZHIuQ+RigKRBC2M/Q2RKqUq\nkVYWKptInwdP6RAgUrVAJEGuQfQ/vdqFAJGqBSIJcg2in/5k7ESASCSPIBIbEEmQaxD95B3R\nDLtIMIkLiCTINYj+qUt+fmIEIlULRBLkuvx92z+7CDEAkaoFIgnCB9Gn9WzwX+MLASJVC0QS\nRAyiT+prB5E2AEQSxPT+Jg2izwa/SDCJCYgkoIu0n59HKv0OWYhUKRBJkOsJWc73I0GkykjP\n9i2LdAjbIzG+H4ma+VSPIFIiEMkk1zkS5/uRIFJlQCSTXFftON+PRE19iFQIiGQSfB/pSL2P\nxPh+JG6RYJIDYgNBJJNcPRtY348EkcpAbSOIZJJLJNb3I0GkMlCbCSKZBA1+EtBFiPX9SOwi\nwSQb5CarWqT5h7JkG46L8/1IEKkI5FaDSCZ0kT5fl79/aQNE8kEWKcAjiGSB3m4QySSmixBl\nyGI+sogEk0zo7QaRTHJ0EVp+m1IIEKkMAQ0HkUzoIl3Jr3WBSC0S0HIQySTgYsPH+KIx6wuP\nzKt24xlVag8hiFSGkKaDSCZR49pZ9jPm2yh8+68Q8ojEY9IbGRnSdBDJJJdIrC8aq1QkRiNX\nJ6jpIJJJrp4N/mt8IdQqEqORFRDUdhDJJJdI1248o+qs724OgCpSmEeJArAqWQFhTXeGSAYB\nIn35RvN2vmgs9V3MVYrEZ2QdBDbd+QyRdEJ6NvguaDteNHb6Sajbi/+I2VpQJEYl6yBYpDNE\n0gjp2eDrx525Z0MGkeLzn1HJSoBIyUT0bLBPTa6Ji+pEYlSyFsJFOkMkFbpIH97RvNWeDdwD\nRNYkEmOoaogQ6QyRFEIeo/CN5r0ZkRhDVUNoyzUhUmGh2nj1ZRaRYtKfMVRFRIl0hkgybyNS\nuEcR2c8YqiYgUjrZ3mp+vw5dGvbX5Nck1SMSZ6yqiBPJ7CoEkdIx3480DsbFMtJqHpECs58x\nVF0EtxtEMskl0rH7GPZF92ty14ZKROKMVRnB7TaKZPS5g0g0vo70UYSYe3+vLRJnrOoIbrdJ\nJL3PHUQiETxk8cC9aZE4Y9VLcLtBJJOQR82f95F+j6RRhPylQ8go0kJUxlBVQ17H0R8h0pkr\n27ckUjfvY0ijCHn3XyHkFMkXljFU3dBXESK5yfeisdeQ+6kDFlNFCk37pfTni1Q9hHXbzWgi\nnZmyfUsinbrnPaH7h3Uns/YN2fDE9+U/W6AWWFyznQxEskMXyX9nKLdIS3kak/rOwFxx2mBh\nvXYaF0WkM0R6EXDVzttXoV2RjMhccRrBv1LCnssgzuuzLNIZIj3J1kXo0/dgegjZRbowRUpc\nz7WwrIllJ/TkPJqkiHSGSAMRIpEuNvgfTA8hv0gXnjiJ67kW+mroFu3EpPNoUlMiFRKKLtLv\nQbhhvobPvCGbfr3uRQGRLixRmNa3NMo6GBbt5KlnxSSzzx1EoiB71HX6oCbOLkLJlBCJB641\nLou8BjZ7JF7eTCaZfe4gEoWuu00fbFO17yfvg+khkESybndXPuSCaYXLojSYz6LLfEN2NMns\nKgSRKOynDg0kkW5734PpIUSL5M2JHPCsb2HUBvO32NSzYbozC5EEua7acY/ZECdSWZMSV3Qd\n1Abzr+DcRWi8M6t3FYJI6dQn0tJxSgYS13QV1Abzr5/oa2fcT4JI5JKKGotvo+AjTaSiJmVr\ng3yoDbawflKnVf1+EkQil2xQpNJHd9naIB9KewWIdNZvzEIkHqoT6ZUXZU3K1gb5UBpsaf2U\nxyj0G7MQiYXM75CNEmn8ByZ5kNsrWCT1xmwPkYi83i+xPGZDeZEsG33Ki6ImJa7qChjt5UMR\naTYJIg3kGrOBj/8ISWpuc+FPSZOyNUIuLO3lQRXpovdwgEgkvrrne8p/7L3oKhQpKEd4yNYI\nubC1lxtNpIvewwEiURBvhTV6rPZVi1Syi0O2VsiE1ERxIik9HC4QiUDomA1cRIm000TCLsmO\n1ESEtdNFOkOkmZDe3w3tkdTEgEgOpAaKEumsdhW6nJOzPaNIeY1603MkiEQisIGsIu1kkc5c\n2f7GIrV01W5niASTbAS2jymS2lUIItF4jVRHG/ubD4JIxhbXEwMi2QlsHotIyqPnr053PURK\nACIt1rJCApvHJpLcw+HV6a6HSAlUJZJxqFLu2C5bM2QhqHVeqWiIJPVwuCg9HCCSHX+nn9pE\nWvwlF9naIQP0xnmVf2biRRdJPHoOkShAJBLZ2iED5MYZy4+ZqIs035gdj/22LhKtg+nvcZ0h\niwNEshyq4NjOhmgc2jpNmaiLdFZFMt+JuSmRup6mw936MsvKRDJyASJZoLWNKC9SUhNpvDEL\nkfrJBIoP1R/arSpSSybNTeNrG6m8lJKmSPPo+hfz5bJbEukFwYcv2ovG2FgWSdvy1sSASBYo\nTSOXl1NSFUnu4nAx34kJkeRJM5+2qay1kvnPtkkVtE1vTQycJJn42kteGXuuns0bs5sW6Q/x\nzafDPOy3zaPqRcJJksliy4zlHLnqvTG7PZGUb/E6QKSZbC3Bzdww9paZy7ly1XtjdtMiJdjQ\ngkjYJan4msu4xmDLVd+N2S2LlCJD9SLh2M7A0y5KOXeu+m7MVitSHqGIJ0iL1COSSxiIpOFu\nF62gJ1dVkS7KjdmtipSmwooiaVngEQkmybiaxSjoy1VVJOXG7EZFShyErn6RsEvSsLeKpWCY\nSOK159sUKRGIJMjWFLxYWsVe0J+rikjqs+cQKZw2RGI3yRoxW1PwYmkue8GFXFVEUp49h0jh\nNCBShl2S3c1sTcGLpVHsBZdyVRFJefYcIllY+3kkokgeW7hF2r0wJ2RrC07EOizVezFXFZGU\n11RAJJNaRaLbwizSbucyKVtbcGJpE0fJ5VxVRJrGBIdIXlZ7sI9FJD6TJodsQbO1BSdmc7lK\nEnJVF2keXR8iOVnrwb50kTh3SWJXZDEpW1twYjaJqyQlV2WRRF8hvasQRJJp9tCOTyTlkK5R\nk4y6O0vSclUSSR5dv16ReIUKF2mtB/tcWzpAFrZjOzUQoXNAhRit5SxJzNWL9GCfEOkMkTRW\nf7CPJJJfFSaRdHMMk7I1Bh9Gg7iLUnNVHkFyJ0bXh0gqqz/YZ9/UQarwiGTugfRfsjUGH0aD\nuIuSc1UVaQeRImhGpHSTbEHa2yXpq+IpSs9Vew8HiBRAGyJx7JKsMra3S9Kbw1M0IFfVHg5G\nVyGItEgJkawbO8wUHpGWf83WGmxotfYVDclVa1chiKTwdVi3ZwOTSIkmOQK0JpJea1/ZoFxV\nugpNPRwgksTn2l2ErFubmOfk6Uu4TGzt2E5rDW/ZsFwVr6toRSQWoULeam575eVEJSItepIo\nknuP1qJIuxwiiddVSF2FIJLA//RsQyIlmOSZu1GRKLUNTWZzWKHzOSjCe4v00d09U1sRKWmX\n5LNQm5atOZhQ28JfNjiZ5WGFIJLB6fjrnlhEJMsGD9YkQST/3gwiiQKSSEqfO4jU53seaQ44\nB1aXwC9SnEm7hTmbOrbTquwvHJHM0/NIu7mrEESayCTS/Fom88NYQ192LmQzIeGpLHnU1rGd\n1hL+wjHJPD2PJPrcQSQSsSLNr2UyP4yEiESRJE6kRY/a2iWpFV4oHJXMQiSpz13dIs0/xLC2\nSPPcJUUKNomgEUTSCsgvTpr63EGkRfKINLx05j9fdvpzmZTwpBkiDE1qkNyoDeEolJjMZ73P\nHTnCe4uUcfATcVa0sEdaQaTdjuhRS7skrb6OUqnJbIwqBJEG3kkkskkBGjUn0uKRHZ9IU1ch\niCT4PbI/2DdfqwsWaSGVaQnvLUjXqKX+dmozuEolJ7P5BjKIJLhzP2reqX/rESnIIkvg+CbJ\njVpbV6n0ZBZdhSCSCfcNWfmfNJGIaU8/5Qm8LNGKSFptXcUYklnrKnQ5kyJsQyTmUYQ6+cPS\nDVl9o/sTmZrwrjLBl8lbObZTW8FZjFWk6Q1kECnPKELS5QtCFyE+kSg9icJv3LYk0vJdJI5k\nFl2FxMtlFyOsLFKUUBEi/VinBi6XjkekhTwmJ7y1QEwHCH2mbI2ShlpXZzGWZFa6CkEkwe14\ntP3cnEhLgw1Fd8iTv2drlCS0xnKWYxZpWJzSVWjbIvX37mr5tS2R/CWjPWrk2E6tqrscTzLL\nXYUuSg+HjYu02tjf5mYPscMo6n1EL/aZpSaO7dSqussxJbPoKrSDSBJfq4qkbPeFLF7IePfg\nC3zP0GZrlRTUxnKX40pmo4fDtkUSV+0KH9q5/x9vZnFQypNUCKSBXZJWUXdBtmQ2RhWCSH/s\nbR6tIhLVDVfG20onetSeSJ6CfMms93DYtEh+mhPJXjrVoxaO7dTV9xRkTGZtVCHjLc0QaaJJ\nkfTiKZcZXJXI1i7RqPX0FORMZq2HA0RyUUwkacsvpHBoygf3USVFzdYusWj19JRkTeazeJXf\nxXxL85ZE6lT0qYHLpfOfMzMXUpiQ8tIcPBo1cGynNpavZAaRRA8HiFSNSJYUDs/5Me93XB7V\nv0tSa+kryZvM4g1kl0ZEChIq4NDuYz8MEPmzt73UvEmRZpMYNbJUI1vLRKHV0leUOZm1V/lt\nVqRr9+/577/i95HsIi1nMDHtWTWyVSNb08Sg1dJXlDuZL8qwQpsVaT6aK92zwS4SJYPpqc+m\nkS1atqaJQdRysWrsyawOK7RVkfbzHon1wb5FdJFeG9+awHw2pFD3Lkmto7cofzKrr/LbqEif\n3X54Eulnv/Kh3XPjk/J3JSwVydY4wWh1tJZxJRaTSMKkbYrUn8brdaWfR7KIRMzfdbAdKWZr\nnVBEHd31yidSr3QV2qpI/c+g0unbOi2rSPpJkit/eYWIph2R7GUyivSQuwqdzT532xDJR0mR\nyOm7EraaZGueQNQq2svkFOkhdRU6m12FIFJBkQLSdx2sNVlcT3LBFLQa2gtlFekxdRU6yy+X\n3ZZI389DO+vYJ+uLtI5HveP6YeAuyRWbHVFDzyLKiLTbrEjH8WJD+Z4N1Yr0qiG1Lvb1Iy6E\nB62C9kJ5RXpMfe52Z7Or0BZE+houf3fd/dh9WaZuUqS5juS6aOsWupxk1Ao6CmUW6dXDYae8\nE3NLIh2GG7Jd19+7g2XqBkVSaqnXhbkyTG0pVdATNbdIvbvP3RZEevYMGv6U7yJUpUhaNXNX\nhqUtteo5SmUXqVde5bdVke7luwgRsnJtj7KLxKKSVj1HqTIiWbsK1SkSRSi6SB/DudEa50jU\nMe/5U9eNpZ5abXJUJ7Upper54uUXqVdf5bctkW777vV0X/Ehi9sTKVd1EptSq52rWAGRHur4\nXJsSqb+dBpHK30eqUCRrRUtUJ60pReVWF+mhvspvUyJ52ZRI9opq1clUn4SWlCrnDVVEpIcy\nPtc5ZE6IFEUzIhXZJaWYpNfNVa6cSKKHA+1VfukFKhAp41vNvZBEqsGjQrukaJP0qjkLlhHp\nIb+B7Ex7lV96AYjkp6RIzpoWq1FcQ+o1cxYsJNLjrNyYpbw4Kb1ABSKN/B5Pd8vP2xHJU1Wt\nRnWZJNfMH6SUSI9JpPF1L8uvqUhfdj0ile8iRNHE9eQCY/JKUV0sV6lELRYqt3OJtJhY7CJp\nN2aXxzJOX3ZFIq3worF4kXp2l7x11apUk0lyvewR1hJJenHS0siR6cuuSaSr7TmKdUVaeJSO\nL30X8ne5UsWq4qya88huBZH681m9n7Qw4F36smsSyUrVIvGZtFRZrVIVmSRXyz73GiIZ73vx\nj9OVvmyI5GXxkdRSqavXKvMlEHormpUyiqwikvS+l90sknNUlPRls9TaBkTizdvFajFDbUWz\nTkaRdUR63Zh9XbzbzSI5RkVJXzZE8mIpoEUplbR6tSoxSamTY8ZVRZpMmkSyP4OevmyI5GNx\nh9RzmESq7XLF2AmqlvvIbi2Rxhuzk0qzSLZHZ9OXDZF8LO+QGEQiVlevWBUmKTVyzbSaSI/p\niG5WybzoAJGS+e+yaBJFpFSTqNVdrFkO6JWqUqR+OqLb7aarDvpFB4iUTJRIljh5U9WxnDK7\npKX62epjllpRpOkyuHh1on7RASIlwyVSkkkB9V2qWi5INfJ4tK5Ij/mG7G4+wDOfr0hfNkTy\n5ZAx2R4pQ44uLqXYLslZTb06zrLrijS9rfksq6Q/X5G+bIjkSSDaDilapNAKG5Vb9SF4ozaO\ncquLNL1kVlZJf74ifdkQyZM8xB1SnEnhFTZrt95IYc66WOq9ukjjS2ZFV4f5Bi3fsrcs0nT+\nKSUD5fs0v/Q9aP4p3TzxvN9D61vi+/CPs/5jmojprx+k8sr3h/m9V74/zO+98v1hfhfnSCPi\n+5TM8veH+b1Xvj8867e4vv7119pvUyL1gYlnzB/yPbS+VrYmUi/E0ZWqTaSH2n4NiDRlgBN9\noieYO4iF2BrbKuiq/o5CUK09iFhydRePZBYnsBaYezYoN2gvTKOjsNZaSdPYbNFYUSRjmidY\nSNrFV9mawZYVIFnEKJMII9eWnjf5UlIuIESa7tDOT/ylLxsi8YgUYFJCla111HUgWcIrUxsi\nzZ3vRLfwy+sHiJQCr0hkk5LqbK1k/J6GyST7kV2FIok+Q/PoKBfl6h1EiiBYJH84YtKl1dmZ\nx/HHawwq2XdIVYo09xmaT5WYhhmCSFwikUxKrbOvrvG7l1SVmhJJPleaTbqYj85GhWaqtQxE\nspBe6dhc95NmkuPIrlqRHtrzFWeO8bogkjN/wk6R+jJvP45N9iVSVHLskGoWSRzivUy6mO9B\nTwidXGuZRkTymBS4Q+qXspyn1nG5TiDapF2bIk29WcXVO8bQibWWgUgaTLWOSXUasTsll0fV\nizT2ZhVX7zhDQyRLgljTxIYn3bhqHZ7odKJMcu6QGhBp7M36XIXn/SSIFMGCSEZSEUI6s42x\n2qGJHkLMTkmaQaspPW/ypeRStovr4IkjSEIkj0jKd1pQeyzWeofleSDBJrl3SE2IpI8SXodI\nc4F+syJZ05y54iF5HkzgTmnn3iG1IVIv+t6po6IwhE4v0G9XJDPN2StOzvI4glSSS+r1JOfN\nqiI95Mvg8UOxQiRHzsScIr3QAmWoOVWJSAI6SPh2SK2IJF0Gl1/fzBI6tUD/FiIp3wMCx84X\nuYQMkFWy7pCC82ZlkfTL4AFdRC8AAAunSURBVKyhx3/lRzgCBl7p2xHJZVKCSH3cXFELyAVN\nJfsOKTjbVxfpdRl8N10GZw3dTy9rUkSS+vf5IvTNixR/ZFeENEtI7Agu2c+QgrN9fZE8l8FZ\nHLWLdBlvXb25SOoP2aoRR5okNHZLLu3eR6RxtzGbFDymsX3C7ItTJP8jHP3biZStFpEweELA\nq9LO4VGjIj0Uk4LHNLZNUHxxi+TpMNtDpNzwmLKMyyX9R6lqjYr0kG8oBY9pbEy4SJ4siuTq\n59e3LlLlp0gDfK4sYXHJlEuqWasiPeYbSrvwMY3VCYYniyKdrb0q+jcQSf0hWy2iYXVlgZ0V\nV/s0K9JrdJT5MvgruYNDuzwRl29cBcx7WP27iZStEilkccaFV6O3EekhXlsxZzslwsIOx/a/\nIccuq12RrCY1IVJhlbzjQsi1almkXntZ5sV7jZp05GbfnYtdkxHhjURq4BTpBSX/k2YmosRt\nWiRxGXweipV0puOYYHqjmuWI8EYiEZNxdYjJHRdiiyI95GxPE0nxRSug7NwtEd5TpGx1YCBN\nIk+MWI9aF2kaJNx//LUkkmVPpP+remYWeLUmRGoKiCRNOEtDsXqOvxwiqedB/jm1Yzy9wNCa\nTYmkmdTMKRInbB61L9I8SLjtDMfuixXSvkxWzihwaUkkc5dk/JCtDlUBkaQJAb7YJZoEWhLp\nrJ9LKSJdLu8kUrYq1EaaR9HZXqNIvcj2QIU891sdIl2UfV/TIini6N+3I1KUSvO87yXSdNHB\n3G1YfQnoSucwzRa6LZH0PZCxh8pWhRqBSPO/IfsVj0gBobVDvAZF0t7UtWGRQlUS872dSI+A\nPtz2CRGhhUrtiXTxi5StBrUCkeQCsSJFh5YuVzQo0s7+ZZMi9QEuSfO8qUh5Q1tVnM6VmhPp\noonkzJUtEeoRRIoNbdmnzdcd3kikbBWoHohULrRxcDip1JxIO/ERO6SZAI8gUmJo+7lSWyIJ\nfYwzpG2L1Ae8qyY6sSDSNEG/jtGkSDt45ICkEURiCa2K9JeQjYk07pJMjyDSi0WNIBJX6Ga7\nCIldksUjiCRYaJXoxIJIeoGGRXKNdp1t8W9Det5AJLNA6yIRjmCASnreQCRrgZZFMjyCSIuk\n5w1EchRoUiTH2MXZFv82pOcNRHIV6JsQSb8SBY+iSM8biOQq0DfxqHloFxhgJT1vIJKrQA+R\ntkN63kAkV4G+DZGWTcq28DciPW8gkqtAD5G2Q3reQCRXgR4ivT98eQORXAX6NxEp27LfAb68\ngUiuAn0jIgU9bgNU+PIGIrkK9BCpdmzbLDFCdN5AJFeB/j1Eyrbo/CRtu7hFQCT+0H0rIvlN\nyrbofLBsO8IiWPMGIrkK9BCJBVvDEuawzBm+ccmVSs8biOQq0L+FSNmWTMbWsIQ5LHOGb1xy\npdLzBiK5CvQQKQlXA/tnyJuSbYZ+MM25Uui+GZHoY+SUhNLA5gzNZjtEchXo30GkbAteZvW0\nsBRoM/SDac6VQvcQKYnV08JSoM3QD6Y5VwrdtyOS06Rsy3VTxbZzFWgz9INpzpVC9xAphiq2\nnatAm6EfTHOuFLpvX6Rsi/VQxbZzFWgz9INpzpVC9xAphiq2natAm6EfTHOuFLpvSCS7SdmW\naqOqbecq0GboB9OcK4XuIVIIVW07V4HlOf83kCc0U4HmQveti5RtoQpVbjtXAbpI/xs/MIZm\nKtBc6L4lkWwmZVuoQpXbzlXAOeF/ikA9RGIM3TcuUrZlvqh627kKQKQVQvcQyUfV285VIFwk\nXSiIFFygb1ukTAtsY9u5CkCkFUL3TYlkqJRpgW1sO1cBiLRC6L4xkfqsHrW17VwFjAm6L06R\n5gLk0A22GER6kUOjNredqwCDSK4bTRDJVUBP0wQKiTSpxLmANredqwCjSORjvoZaDCLNhGr0\nptvOVSCDSIvHfA21GEQKpYoGzh3aUmCeEOzJskiuY76GWgwiUamqgXOH9q0vRCoXOjxNnUCk\nqND6aUj6skuIxF9reotVGTo8TZ2sLlKVDeyc4MrV9GXHexIuEl+tV90Y6QXoaboIp0hdJ0fb\nlEhzzkYvew2R0mu96sZIL6CmqZq/gTCK1KnhXjWU66xQdQPrP6QnM0O2FwjdxMZgDN3LImn5\nGwifSJ0WTxWprQbWf6gq2wuErnpjMIbuJZH0/I1Kfw7sIiVumsUCvjkTU5I1mTNk+5qhnRvD\nOWfG7RxfoK9dpP8AaAcjf6PTn4PEikhwXUhpLWQTlWwjZEREiPQ2IZuoZBshIdJAE5uqjq2P\nkGwRIdLbhGyikm2EhEgArEM1IiXe0AJgXWq5IZvYxQKAlamlixAA2wUiAcAARAKAgfVFmo9M\nO/FJ+5kp5Ly2fCHnb+EhLSvIud5SSK71VmsZu97yLEwr7oiYst7BrC7SfK1EXDTp7D+nhhwa\nVJ2+Zi23ut7yLMr/kuJDOiKmrHc4a4s0r6q04Tv7z4khhz+dMn3NWqqzmoGSK6nUlGe9eWs5\n14utKdWIKesdQTUiSV/5NpUc0pZiqSHja6nOwZ2iyje29TZboWqREtY7ggpF6pSfuz60ko6Q\n6iS2kJG1VCLq1WIQyaIAb8j0WnbWFU9pSqNGsds7grpEkvce89fQM0VnDomUYAw5TQsMaVlB\nRpHkkFqNuULqmyk8pEOkhKZ0iLSJiw2LIqXuPiwisYbU65wWUVtvjpDqJLaQUp3jQpqrm7pH\nMtXe6h6pU/+k/5/ZufFXD2n+710PNPxvNCHrReOJH2sKqW3rhA3uiGiGzkpNIslZNH1MzHoz\nE5hDWhcTGLG3ipRUSSmkFocpZGeUCgzZjTCKpEacJ21OJPN/+MkiySHN2AwhO2N6aET1A8t6\nax96hvXWQ0avtzoL6x5Ji2iGzsraIj0r4MjQTvmQHLKX25QnZGdMD47YjR8Y11sL2TOst6WW\nqSG1zwkr7ojYJ6x3OKuLNF1SmfbLnfpzH3PJxRFSbDSmkPOxRExIbQX19e7iK6mH7HvmkGK9\nk0JqCrB1EeLY3sGsLxIAbwBEAoABiAQAAxAJAAYgEgAMQCQAGIBIADAAkQBgACIBwABEAoAB\niAQAAxAJAAYgEgAMQCQAGIBIADAAkQBgACIBwABEAoABiAQAAxAJAAYgEgAMQCQAGIBIADAA\nkQBgACIBwABEAoABiAQAAxAJAAYgEgAMQCQAGIBIZXi9COVwvesTfhZnvX903XWK4ppiDVvg\nbSZgBE1dhum9jPub+vtheQOc/mb7nKK4pmgcOktxkBE0dRleOX07dkfb7wuz3iKmBNQNMIAG\nL8OU2Yfux/o7YVaeKSAPaPAyTJn9030Mf/8OyvbXfjziG37/OnT7L3mGvx8OX1OJTorytxs6\ndftPaYo073XfHW9z2Dn0M5KYE/ADkcowuXDvDn3/+XLgKkQ6PT9Ih33H6QeLSPvueW40TZHm\nfc61vysizZHEnIAfiFQGzYXvvv/u5teF/+2njvf+fhSHfd/d/l//bz+Ukw/TXjP/lf0afHxN\nkeb9Hj5+vASdikuRpDkBNxCpDIpI0qfX11M3XBa/d6dp2unp1M+4I1GidN2vMrM072mYdO/2\nskhSJGlOwA1atQyaSLefz6Mk0nRxvNOKS0Xm78KR6c88rz6/9Js2J+AGrVqGKX1vz53McdYG\nIr0JaNUyTOn7PZzBfHSHr5+bIpK9OE0kYyEQqTxo1TKI+0i/4xdZpJN2d2k+szn1iyJJ8x49\n50gniJQVtGoZlJ4Nw1n/v+kcaeic8Ly01n+Jiw3eq3ZzvOcfad6v4bLc9XXV7tZbrtr1PUTK\nBFq1DEpfu+v45XfYQw07kPGkSeqIJ+7+LIkkzzvdRxrDGveRtHCAD7RqGV7mHMeboR9/H3+f\nh1u/h6dIQ/eD7kPuOfe1H/sjLIokz/un6Gn49Ao79mzYzz0btHCAD7QqAAxAJAAYgEgAMACR\nAGAAIgHAAEQCgAGIBAADEAkABiASAAxAJAAYgEgAMACRAGDg/8HzC52+plHoAAAAAElFTkSu\nQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 300,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "options(repr.plot.width=7,repr.plot.height=5)\n",
    "\n",
    "scaling_parameter=max(result_withcut$total)/max(result_withcut$upper[!is.na(result_withcut$upper)])\n",
    "range = c(as.Date(\"2020-02-01\"), (datestar-17))\n",
    "adj = 0.8\n",
    "options(warn=-1)\n",
    "\n",
    "figure_data %>% \n",
    "    ggplot() + \n",
    "    geom_bar(aes(x=onset, y=value, fill=subject, group=subject),stat='identity', width=0.7) +\n",
    "    scale_fill_manual(values=c(\"#FAAB18\",\"gray35\")) +\n",
    "    geom_line(data=result[!is.na(result$Rt),],aes(x=onset,y=Rt*scaling_parameter*adj),color=\"#1380A1\",size=1) +\n",
    "    geom_ribbon(data=result,aes(ymax=result$upper*scaling_parameter*adj, ymin=result$lower*scaling_parameter*adj, x=onset), \n",
    "                fill=\"#1380A1\", alpha = 0.4) +\n",
    "    ggtitle(\"Entire Japan\") +\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",
    "    scale_x_date(date_labels=\"%m/%d\",date_breaks=\"10 day\", limits=range, expand=c(0, 0)) +\n",
    "    scale_y_continuous(limit=c(0,640), expand = c(0, 0),\n",
    "                       sec.axis = sec_axis(~./(scaling_parameter*adj), breaks=c(0,2,4,6,8,10), name=\"Effective reproduction number\\n\")) +\n",
    "    geom_hline(yintercept=1*scaling_parameter*adj, linetype=\"dashed\", color = \"#1380A1\", size =0.7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating Rt excluding the cases with unknown illness onset date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "entire_Japan %>% rename(domestic = domestic_org, imported = imported_org) -> df_onset_wo\n",
    "\n",
    "## adding extra 10 days for the stability of back-projection procedure\n",
    "temp_lastdays <- matrix(NA, ncol=1, nrow=10)\n",
    "temp_lastdays %<>% as.data.frame() %>% mutate(onset=as.Date((max(df_onset_wo$onset)+1):(max(df_onset_wo$onset)+10), origin=\"1970-01-01\"),\n",
    "                                              domestic=0, imported=0, total=0) %>% dplyr::select(-V1)\n",
    "rbind(df_onset_wo, temp_lastdays) -> df_onset_wo\n",
    "df_onset_wo %<>% mutate(time_onset = 1:nrow(df_onset_wo))\n",
    "\n",
    "\n",
    "## incubation period\n",
    "K = nrow(df_onset_wo)\n",
    "incubation_probability = plnorm(1:K, inc_fit$meanlog, inc_fit$sdlog) - plnorm(1:K-1, inc_fit$meanlog, inc_fit$sdlog)\n",
    "inc_pmf = c(0,incubation_probability[1:21])\n",
    "\n",
    "\n",
    "## back-projecton of imported cases\n",
    "sts = new(\"sts\", epoch=df_onset_wo$time_onset, observed=df_onset_wo$imported)\n",
    "bpnp.control = list(k = 2, eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                    Tmark = nrow(sts), B = -1, alpha = 0.01, verbose = FALSE, lambda0 = NULL, eq3a.method = c(\"R\",\"C\"))\n",
    "sts_bp = backprojNP(sts, incu.pmf=inc_pmf, control=modifyList(bpnp.control,list(eq3a.method=\"C\")))\n",
    "df_onset_wo$imported_backproj = upperbound(sts_bp)\n",
    "\n",
    "\n",
    "## back-projecton of domestic cases\n",
    "sts = new(\"sts\", epoch=df_onset_wo$time_onset, observed=df_onset_wo$domestic)\n",
    "bpnp.control = list(k = 2, eps = rep(1e-4,2), iter.max=rep(1000,2), \n",
    "                    Tmark = nrow(sts), 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_onset_wo$domestic_backproj = upperbound(sts_bp)\n",
    "\n",
    "\n",
    "## normalizing the back-projected cases\n",
    "df_onset_wo$imported_backproj[df_onset_wo$imported_backproj<=0.01] <- 0\n",
    "df_onset_wo$domestic_backproj[df_onset_wo$domestic_backproj<=0.01] <- 0\n",
    "\n",
    "df_onset_wo %>% mutate(imported_normal = imported_backproj/sum(imported_backproj)*sum(imported),\n",
    "                       domestic_normal = domestic_backproj/sum(domestic_backproj)*sum(domestic),\n",
    "                       total = imported_normal+domestic_normal,\n",
    "                       time_onset=0:(nrow(df_onset_wo)-1)) %>%\n",
    "dplyr::select(time_onset, onset, imported_normal, domestic_normal, total) %>% \n",
    "rename(t = time_onset, imported=imported_normal, domestic = domestic_normal) %>% filter(onset <= (datestar-1))-> dt.backproj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "today = as.numeric(datestar-as.Date(min(dt.backproj$onset)))\n",
    "\n",
    "## start estimation from 1 Feb 2020, as there is a huge uncertainty in the early phage of epidemic\n",
    "Start.T = 29\n",
    "\n",
    "\n",
    "## reporting delay distribution from infection to lab confirmation\n",
    "infectiontoreport <- function(t, tau){sum(convolve(onsettolabconf(1:1000),\n",
    "                                                   rev(incubation(1:1000)),type = c(\"open\"))[1:(today-t+tau)])}\n",
    "\n",
    "## pre-calculation for reporting delay to reduce the calculation time\n",
    "precal <- matrix(0, nrow=(max(dt.backproj$t)-Start.T+2+1), ncol=(max(dt.backproj$t)-1))\n",
    "for (m in (Start.T-2):max(dt.backproj$t)){\n",
    "    for (n in 1:(m-1){precal[m-Start.T+2+1,n] <- infectiontoreport(m,n)}}\n",
    "\n",
    "delay_precalculation <- function(t){sum(convolve(onsettolabconf(1:1000),\n",
    "                                                 rev(incubation(1:1000)),type = c(\"open\"))[1:(today-t)])}\n",
    "dt.backproj %<>% rowwise %>%\n",
    "mutate(imported_delay = imported/delay_precalculation(t),\n",
    "       domestic_delay = domestic/delay_precalculation(t)) %>% mutate(total=imported_delay+domestic_delay)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "## estimating Rt\n",
    "est.t <- list()\n",
    "est.CI <- list()\n",
    "\n",
    "for (TT in Start.T:max(dt.backproj$t)){  \n",
    "    \n",
    "    dt.backproj %>% filter(t <= TT) -> dt.backproj.T\n",
    "    \n",
    "    llk <- function(param){\n",
    "\n",
    "        llk <- 0; Cs <- 0  \n",
    "        t=TT; R = param\n",
    "\n",
    "        Css <- rep(0, t)    \n",
    "\n",
    "        for (tau in 1:t-1){Css[tau] = (dt.backproj.T$total[t-tau+1])*generation(tau)/precal[t-Start.T+3,tau]}\n",
    "\n",
    "        Cs = sum(Css)*R \n",
    "        Cs[Cs<=0] <- 1e-5\n",
    "\n",
    "        return(-(-Cs+dt.backproj.T$domestic_delay[t+1]*log(Cs)-lgamma(dt.backproj.T$domestic_delay[t+1]+1)))\n",
    "    }  \n",
    "\n",
    "    param0 = c(0.7)\n",
    "    opt_est <- optim(param0, fn=llk, method=c(\"L-BFGS-B\"), lower=c(0), control=list(maxit=10000))\n",
    "    opt_est$par -> est.t[[TT]]\n",
    "\n",
    "\n",
    "    ## 95% confidence intervals using the profile likelihood\n",
    "    ci_pro <- matrix(NA, ncol=2, nrow=1)\n",
    "\n",
    "    CI <- function(par_CI){return(2*(-llk(par_CI)+opt_est$value))} #llk()\n",
    "\n",
    "    par_CI <- seq(0, 10, by = 0.01)\n",
    "    \n",
    "    logLik <- sapply(par_CI, FUN = CI)\n",
    "    as.data.frame(par_CI) -> par_CI; as.data.frame(logLik) -> logLik\n",
    "\n",
    "\n",
    "    cbind(par_CI, logLik) -> data_CI\n",
    "\n",
    "    data_CI$logLik[data_CI$logLik<(max(data_CI$logLik)-3.84)] <- NA\n",
    "    data_CI %<>% na.omit()\n",
    "\n",
    "    min(data_CI$par_CI) -> ci_pro[1,1]\n",
    "    max(data_CI$par_CI) -> ci_pro[1,2]\n",
    "    \n",
    "    as.data.frame(ci_pro) -> ci_pro\n",
    "    colnames(ci_pro) <- c(\"lower\",\"upper\")    \n",
    "    ci_pro -> est.CI[[TT]]\n",
    "}\n",
    "\n",
    "matrix(unlist(est.t),ncol=1,byrow=T) -> est.t\n",
    "matrix(unlist(est.CI),ncol=2,byrow=T) -> est.CI\n",
    "\n",
    "cbind(est.t, est.CI) -> est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "est %>% as.data.frame() %>% mutate(t = Start.T:max(dt.backproj$t)) -> result\n",
    "merge(dt.backproj, result, by='t', all.x=TRUE) -> result\n",
    "colnames(result) <- c(\"t\",\"onset\",\"imported\",\"domestic\",\"total\",\"imported_delay\",\"domestic_delay\",\"Rt\",\"lower\",\"upper\")\n",
    "\n",
    "result %>% filter(onset<=as.Date(datestar-(17))) -> result_withcut\n",
    "\n",
    "result_withcut %>% dplyr::select(onset, domestic, imported) %>% gather(onset) -> figure_data\n",
    "result_withcut %>% dplyr::select(onset) -> TIME\n",
    "rbind(TIME, TIME) -> TIME\n",
    "cbind(figure_data, TIME) -> figure_data\n",
    "colnames(figure_data) <- c(\"subject\",\"value\",\"onset\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAAJYCAMAAABvmDbGAAAAQlBMVEUAAAATgKEzMzNNTU1Z\nWVloaGh8fHyMjIyampqemk+hzNmnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD6qxj////s\n5NGUAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2diZaqOBRF47OdyypLi///1ZY5\n83gDiZy9VvdTCZck3l1BCMAaAEAybO0KAPAJQCQACIBIABAAkQAgACIBQABEAoAAiAQAARAJ\nAAIgEgAEQCQACMgpEuNRln7vh0KesUirFrOZfhlFRRZqDFiQ9UQaPoJI4BPIK1L80rTS0bhF\nyr0VUCcQyXszEAmYWVKk9/vnibHzbzPt9g2FGHsd2On94nXbs/35aYz1/V6/D9B99H2YS0uL\npi3NzMHP7Kv75MEO7///nvfscOc2M1Z8+Ldd772cWybFf74DHB9ie8UY4gr9Z8/9UROqrePt\n2Vbt3G+c9e07s4ehWVz1hU543Q6MHe9y26UFgIZlRXruO31+VZHeCXBrc6v/+McQ6zis9ug/\nunOlpUXzlib44Pt+yb7N0p/+4+O8GVGCYfldEImP/+hf/1hF4lfoPnvu909lybCt/U9bt279\nW9svbbm9oVl89flOGJt7kNsuLABELCvS4f038ffY/63lDza8l7yaV5s899f7jzLbv7Sx7uz4\nzr3ncUz69u3r3KaYsojf0gAfvB+KTuze/slnX6+27L3Ri7Rn53fl3n/2eZH4+GIBqfHjWnLT\nB4/EJa8hVFvHW/8HYoj6aH3SNYuvvtAJR3Z7N/S94Ca1XVgAiFjqqF3/vtud+BWydUjMLmnu\nfYHma9jzEmK9/3fod3NeY4BOml4HeRG/pR4xeLtz990F+OrT8rdTSyPSfUjbkyASF/+72ymd\nCvAV5iJKTX+PYo9GWfI1hfrq1WmDdwVubXFds/jqC50w7pe2f2eEtvMLABXLivQcPp//P/1G\nal+fhtq8BkWEWMqbYTfunXEnddFTWUsKvmff/Zhw5HeUNCKdhuUPcdeOL/DgC0gV1jf9Z9p9\nFZYcp1DvOu7bRD+wV2d490bXLKH6/LbfPzp/xoFdaDu/AFCx8MEG7l9RpOG18bTT+MnrcT/t\nxXSf/sSqi/g4UvB32rNvpZaa1aXNCJ/pCtgaO72bf+eYQ93eRX7e401r1283POmbJfTT3And\nb6fD7aG0nV8AqKhLpMdBGOD4ZdpFFpHef6UPai2XEml/HAddc6jHez/s9B642l3Hr26ocook\ndMKjP/RweMpt5xYAKgoTyR7rPYocvr5/X2JyDlmnW2T50/1k0xE+pYr5RXo+h987tlD7w7P9\n99WeHNj7NEvshPeaP+2B7rPasdMCQEVRIo0/EXhOXLHD8MNiDDD+eDmpi9QaSMGP7Kv/tW3+\njdSl8fQTiD9GIsS3/kZ66kVqjyL0vwuFJcJvpPcw9N0l+/H9k+pmbtZcfbETmrkCuo59Wv5s\ngWCKEmk8aMUfPuiPbj3ZcS42ndEZD6d9q4vUGojB24Nd/QGvrz5JX8JRu9852HhQ7mwS6Xuo\nx1kSiYuha/qpjyss4Y/adY3uzPhmB3HwFJs1V1/ohEMvaveh0HZ+AaBiTZHm5BiWvPbs9Gwe\nt/HoVMuZHV/tqZE2sw5t0rzei7tUeP/Tn3Q5aBapNRCC/7L+kNij23O6N+J5pPc2f5vXffit\noT2PxMfXnkcSY+j+huy7gx3CEv48Ule/buGLST/nhGZx1Rc64d5OtmhjnaW28wsAFXlFEn/k\ninlw6j8Vlgxn6fmZDcN5+G4v7Ltf+nUYfnr3Uwq6g9jyIn5LjRr80I0XvU7qzIbfoWD/7ilM\nXNC05Kmb2SDG0In0YHtlidgBwzmi5jicV7M166h0wjDNoXNSiMsvAESsJ9LzzPppL9ySbkoY\nOws79F254SR8e7zp9Bj2zt4rvUeqw5d+EbelkTn4ONeuz9Tf9xaEuXbdNufJc91cu0djFKmb\nIHeQ5tqJMbR7tV/8cQAu1P42pPhzEOpbOn8kboirvtAJzf3Ipu4ROpZfAGioeD+5tH38F2av\nbZjCkjGEUkQaDr9Nc4nAFikkGWMoRaRT/+P9zjSTdcBWKCQZYyhFpPGyhH7GEdgmhSRjDKWI\n1Ly+DoztzxiPtkwpyQhA1UAkAAiASAAQAJEAIAAiAUAARAKAAIgEAAEQCQAC8on0D4AKIEr3\njCJliwwAGRAJAAIgEgAEQCQACIBIABAAkQAgACIBQABEAoAAiAQAARAJAALWF4nNtyuUX3RA\nJFABq4s03XB0/G9+0QORQAWsLRIb/6++GIBIoALKEKmBSKBu1hdpfGDJGIYXiXJ+OgAZWV0k\n4ecRRiRQKauLNP4fIoGagUgAEACRACAAIgFAwNoiac7D4oQsqI/VRcIUIfAJrC+SC4gEKgAi\nAUAARAKAAIgEAAEQCQACIBIABEAkAAiASAAQAJEAIAAiAUAARAKAAIgEAAEQCQACIBIABEAk\nAAiASAAQAJEAIAAiAUAARAKAAIgEAAEQCQACIBIABEAkAAiASAAQAJEAIAAiAUAARAKAAIgE\nAAEQCSzIX8valcjCpkS60oUCUUAkF3WIBJNWBiK5qEQkmLQuEMkFRAIeQCQXtYgEk1YFIrmo\nRiSYtCYQyUUNIl2vMGkd/kaBIJILiATMQCRvKhIJJi0ORPKmJpFg0tJAJG8gEjADkbypSiSY\ntDAQyZsKRLpCpLWASN5AJGBGEenv44yCSCA/EMkbiATMQCRvIBIwA5G8qUskmLQsEMkbiATM\nQCRvyhfpCpFWAyJ5A5GAGYjkDUQCZiCSNxAJmIFI3lQmEkxaFIjkDUQCZiCSN8WLdIVI6wGR\nvIFIwAxE8gYiATMQyRuIBMxAJG9qEwkmLQlE8gYiATMQyZvSRZI9gkhLMPoCkbyBSEAFIgUD\nkYAKRAoGIgEViBRMdSLBpAWASMFAJKACkYKBSEAFIgVTuEiqRxBpASBSMBAJqECkYCASUIFI\nwdQnEkzKj1OkzxEKIoF8QKRgIBJQgUjBlC2SziOIlB+IFAxEAioQKRiIBFQgUjAVigSTsgOR\ngoFIQAUiBQORgApECqZokfQeQaTsQKRgIBJQgUjBQCSgApGCqVEkmJQbiBQMRAIqECkYiARU\nIFIwEAmoQCQNTEReSlstDohULxBJA0QCoUAkDYM7/T8QCXgAkTTUK5LJI5iUG4ikASKBUCCS\nBogEQoFIGiASCAUiaYBIIBSIpOETRYJJeYFIGgR3IBLwACIFA5GACkSyoYxG3YfJNTEBkepD\n9gQi6ahNJItHECkPEMmCMNPu9JKX0laLAyLVx9ZFstogzlk9y0vJazaSVSSYlIWNi6TdZ1N5\nF3vWdNQOIi3OtkVinja0EkEkYCFYpPqFihTp6yZ/SFclCYhUH5sWifna8ElH7SBSFrYo0pv+\nlUukjzxqB5GysEWRxhes2eRRO5iUgw2LxKb/OansqJ3dI4iUgy2LpL+niY7KjtpBpOXZsEgd\n/jZUdNQOIi0PRLLyfWoPM/xol0EkMAORbBzHI3a6hdWKBJMysHWRrNzZ/uf92+h1ZHfNUogE\nZiCShQP77Q4yvNhBsxQigRmIZKE7UKc9YtctJauRDESqD4hkYRTpxfa6pWQ1koFI9QGRLJzb\n30b4jQQ8gEgWnnvWn7c96pYWK5LLI5iUAYhk43lqRartPBJEWgGIFA1EAjMQKRqIBGYgkpV7\nO7fh9K1dBpHADESyMU4RwsEG4AAiWbix4+P9z+PI5JnfLRWLBJPIgUgWGOsvL6/shCxEWgGI\nZGGaGVTXFCGItAIQycKJ3doh6XXW/kgqVSQPjyASORDJwnM/zLbbPzVLIRKYgUg2Xrf2x9H+\nJt+Jq6NmkWASNRApGogEZiCSD1UdbIBIawCRLDwO3P0hlYtkIRKYiRepWqNCLjXnkeeAQyQw\nA5EsMPYcX+iW0lRHA0SqD4hkYT9OaPg8kWASMRApmkJF8vMIIhEDkaKBSGAGIlngjjTUcxN9\niLQKEMkCRAK+QCQvdHOEIBKYgUgefB8rOmrnKRJMogUiuXicddMaGogEeCCSledt/7boS3cV\nRZJIw7rT8wLFBwdCpHrw88RZ4LNF+j4YLWqSRBq0mZ4GLT0WGiLVA0Ty4D1O6G76PS6NrgHj\nHwTNuBcDi4gEk0iASB78HBnb3+hHpJwieXsEkUiASF68vsx7d9EisQYifQwQyZfHjfqonVWk\nfy2xkRuItDQQKYAf0vNIwiEGjEiVA5GCIJzZwMT/rScSTKIAIiUTK9I8e28MA5HqBSJ58X1i\nWR40hhHpU4BIPoxPozjpFqbPbMhxQhYiLQtE8uDO9j95Hsacb4oQRFoWOpGqEyrkLkK/3XSe\nVz2TVgM8gkkUQCQPhjt/13TzE4i0MBDJg1Gkip6PBJEWBiJ5cG5/G2X5jWQHItUDRPLguWf9\neZ96niELkRYGIvnwPLUiZTiPZGcpkWBSOhApGYgEIBIBEAlAJAIgEoBIBJQoUphHMCkdiJQM\nRAIQiQCIBCASARAJQCQCIBKASF4IT6NQl9JVSWI5kWBSKhDJA4gEXNCLVI1QwWn6OJ5qeawL\nRFoaiBRAPRf2QaSlgUghYNcOGIBIIdx0dz/5CJFgUiIQKZkCRQr3CCIlApG8yHVfOzsQqR4g\nkg+B97X7ObU/pk6m58B4A5HqASJ5EHhfu+NwJ+J9qkkQqR4gkgdh97W7s+OrFenOzkn1W1Yk\nmJQGRPIg7L52e/aa10gCItUDRPIg7L52k3HpIsVmN0RaHIjkQdh97Q7DiPSr3REMYVGRYFIS\nEMmDsPvaDb+Rfva2R6F7AZHqASL5EHZfu9NwsFyrXQgQqR4gUjL680js9J0cOVakKI8gUhIZ\nRSpdqAqmCC0qEkxKASI5YCK6ErTV4oBI9QCRHESI9Lq1R8n3N91VgEFApHqASD4c94/3/3/2\nX7qFskjP/XAWKX2K0LIiwaQEIJIHp3aKUNOeGbpplsoiHdm5HYteN/0c1wAgUj1AJA+mPTqv\nXTt76RAgUj1AJA8O04jkM0WonWvX8oJIGwIiedBdRtH+RvLatbuxY/uL6nHUlg5haZFgUjCB\nnjgLfLZI9rkK+uuR1pzZAJEWAyKFnZD9OZvnKqg7cN2F6cfUmXbRIkV7BJGCgUgfPLMBIi0H\nRKpDpKjUjhcJJoUCkaq49zdEKh2IlFGkr4OldAgQqXggUr6b6H9ZtQtheZFgUiAQKd9N9NOv\njB2BSMWzpEjTB4WR6yb6yQPRBEQqHoiU7yb6J5Z8/cQARCoeiJTv8Pdz300RImAFkWBSGBAp\n5ib6fjMb7Mf4QogTKckjiBQGRIq6ib7XXLuSRdrtIBIpEClu9rfXTfTJoBZp1wKRKIFIISLt\np+uRFn6GLKlIuwGYRAhEyneFLOHzkQhF2u0gUgYgUtwVsl4jEt3zkehEGh1yiwSTQoBI+X4j\nET4fiUykaSjCkEQLRMp31I7w+UiUIimvIBIBECniPNLR9zwS3fORYtI6VSSYFMAaIpUmVK6Z\nDZTPRyISibMH+3akQKR8IlE+Hykiq+0DEoYkWiBS4M1PAqYIUT4fCSIVDkTKeDsuwucj5REJ\n+3ZkQKQQkb76w98PvxtE0pFFJAxJlECkuClCPrcspoNKJFEdiEQIRMozRcj9NKUQsomEfTsq\nIFKISDfvx7rUIBKGJEIgUtDBhvPwoDHtA4/Uo3bDL6rUGUIQqXwgUuR97TTjjPo0Ctv4FUI+\nkWASEWuKVIpQuUSifNBYFpEwJNEBkfLNbLAf4wshRiQPjyASHRApn0g3NvyiYtpnNweQUSSY\nRANEChLpbrubt/FBY6nPYs4mEoYkMiBS2MwG2wFtw4PGTj8JdeuBSMUDkcJmNtjmceed2RCa\n0b4iwSQSIFLUzAb90uSamKARSSeNh0gwyQeIFCLS2Xo3b3FmA/ENIlcUCSZ5AJHCLqOw3c27\nUpEwJFEAkSp59GUmkTAkEVGESCsLBZFgUjIQKeNTzV+3dkrD/pb8mKQIkXyV8dq3g0hOIFLO\n5yMNN+OiuNNqLpEwJNEAkfKJdGTndix63ZKnNqwvEkwyEJ/t2xbpfvS/ixDt7O+MImFIigci\nzeS8ZXHLq2iRMCQlAZFmQi41784jPY5edxGylw6BQiSjL54ibcSk0HZCpJmQw9/jGON1FyHr\n+BVCbpGihqSPVCy4ORBpJt+Dxvpb7qfesDizSN5DklAFm2LVEtEeiDTjL9KJdeeEXmftIFPU\nCdksIo21sAhWMzHtgUgz/iLZzwxVLFKASUa7qieqQRBpJuConXWuQmaRgtJVq4vJg2SRPsGk\nyPZApJlsU4S+bBemh0Agkk0WiNTIPea9mprMl8tlbZGmD5YlQiSvgw32C9NDyC8STIpskJzM\n10vLFSJZeRxmN9TH8KknZNOP1/VkFglDUhN9AEVK5usg0vUKkSzwHjEm39TEOEUomfJFqt6k\n2BYJydwKNIr03sGDSCYYe44vdEul9yfrhekhLCDS1ockApFGgaZ/IZKJ/TihwUuk5952YXoI\n6SI5TMGQFNsim0gXiBSFumtXzsEGiOQgtklzMl9UkS4QKYa6Rdq2SdFNmpL5ApEi79ngfBoF\nHcEi6UyxWgCRIts05uxFK9IFImmBSDayNX8Bots05OzVINIVIoVTuUhbHpLiG9Xn7BUiNXlE\nYiKJkZNFcmuy7SEpvlFdyl7MIiknZiFSR/98Cfc9G2oUacNDUnyrOo9sIl0gkoagezbQsYBI\n2z7cEN+qv+H8EUQKEenOuueU/+hn0W1epGpNSmjXeNjbLNIlSofPFml+KqwyY7XJL1JAnkZY\nsuV9u+im8b4YRbrE6EAi0sJC5btnAxWhIukscaYKRIoBIs2EzP6uY0SSv20vRzY8czWhyapI\n726URLpAJIl1fyNlFmm7Q1JKixWRdj0QycqqR+0gUiZSWiyLtJvhFkAkmf5OdX73/qZjMZG2\nuW+X0mJJpE4gTiWIFEHtIm31ivOkFosiTQPRaJI6eRUiOfkEkbZoUlKDBZFGj8aDDm13QiQN\n9kk/xYrk7QdE8kb+bTR6xB/+5g86QCSeSkSSv3SIZCWqoapI/Ah0ncyCSBYexzVuWbyMSNu7\neXFcQ/UiXUWRLuOuHjcLHCJxvLQPs6xfpE0OSXHtVEQaDtIpIu0gkoUP3bXb5JAU107d+SPN\nFCHufBJEUrn7PWiMjMVESh+SsvVBLiI7QHf+SDfXDiLpmI81fOmW0lVJAiLlQ2m/Xw+IIs0D\njzJpdTRsvJwCIjWzSAedR8WKFCLH9vbtlOb7dYEgknD+SBVpB5GCyC6Sb5Kqcvh7sGmROof8\nTOJFmjzSXkYxTxXakEhJ91coRaQkN5KHpGydkAm17V6DEieSMNlbcz3SNFVoNZEWEmoWiTUp\nOnyESMlDUrZOyIPY8t38ytELgie7q0Ok3bZEYtz/I/gYkbZkkr7hbpP4n0Dz3NSrTqTLZNpW\nROrZuEgbG5L4ZvPtdqk0H63jrpYwi7TbhEhvuLefJlKoGKlDUrZeyIKx1Q6TZo3Go9tmkeY5\ndx8ukvDu434jBXuxVZE0/WB7EPyskSjSnMzqYT2I5MMHibQZk6y9pO2InYB6zwY+mfkTTX2B\nZjMiJdjwKSJtaUji9dB2hGZ/jxuN1Dl3YjIrZ2ybrYhkl2Hl65H8UjRdi8QhKVs3ZMDRSUJP\nWEei4V85meU5RM1GRHK4sBGREoekbN2QAVeLdyr9Ar1IjZLM8uTVbYjkqcJaF/ZFiRQjxXaG\nJGeDtRZdDSJpklmeHr4JkbyfwbLShX2LibSZIcmjvTqLrnqRtMk8FhiPk29BJH8q2rWLFSnB\npGz9QE58J2lEMiTzLNIOIkmsdGHfciKlDUnZ+oGc+Ob6i/Q3FIBIMytf2LesSFswKb61qkjm\nZJ5E2sknmrYu0koX9vllqKpEjAtbEym0iYpIlmS+8JdTQCQnZYhEZETSkJStI4iJ6KPBG0Uk\nazKL9+uCSA6KFCnahxSTsnUEMRFNNYjkSGbucgqI5OSzRNqCSRF9lCjSDiJ13A+rzmxYVCT/\nW++rn2XrCVoi+kgvkjOZ58sprsLhvYVFymuUv0hfmaYITQGnwOIW1hHJa13lTKVvPUsgoo+0\nInkkM0Ti0T/yciTpikCmf9Gznkj2lbXn/D3rWQC+zeTRieSTzNN1SVdhCsRGRbJPH4oVabpV\nhPpiYCWRHCmmnYPmX9H1iekijUh+yXwZhySI1DRn9rIsTfuNlEektElzNpM4hbYtkmcyjyLt\nINKb0/FhXliESGq6+yeJBqNJwkCkFkrqiqVwNFGLLJJ/Mo8mDT+qNi1SxuuR5l9FkkjtrVnW\nE8mUZuL+nFoopSsWI6aLJJFCkpkfkrQFINKwNKUSRpG6Gvrnp5rwAVmiQf8TSP5MKZPSFUsR\n1UWCSGHJ3K0JkdykiDQdqytNJO0Bbu0n4gcJfbEUUV3EiRSczJf5Aj99AYjUkXzflDJFUgcl\nwxglvI/vi8UwN8bCKFJMMl/mC/z0BSBSR/KNvOhFovDo6v5FpPkwui+WI6qLhokMcck8X+Dn\njPDZImX6jcT4F4knZJX0DskSIzsJfRH+bWxfLIi+ixp1Aj2P37wEQ4FuLOtEukIkepG4iK4p\nQu78tGZ3Ai6P6jNpbpem2sZ+SEtmiCTwOK51YZ8zPe3JnYJDo4pFMtRa18bUZB5E2kGkntdK\nl5qvKZITybFsnUGFrod0JcSlicl8HW7ecFXPQG1RpLXuIlS0SNacLBBNpQ1l+IU0Iu0gUsda\ndxGqSaTiTRrrvLPWWF6SmszDNecbF2ntuwiFirSoR5Xt26k95LUakUg715movCLlESpCpB/t\nUrIayfzzTU81tbNpo6FukfxWS07m6eYNWxZp4Hk86j6GSDWL5LkaRDITc9TupvkUItVkktxD\n9tKEyTxdKQuRVjxqV7RINQ1JU4XXEOkCkXruEElLjSL51ZYymSESd9SuzF07a14vgLhvl607\nKJB7yFGcXqSdfJH6FkXa6zyCSK7TmwWh1NdRnjSZx0vOtyuSnUVEsn7jtrRehGr27abqQiQ6\nNimS/WKBWKrZt5M7yFWeNpm5fTuIpFKTSJrSJNQyJMm1dZWnFum6bZGYiLyUtlocniLZs5pH\nX56ASkRSautagTiZIVJNIrk8ymFSJedklR5yrZBDpJ14s+MtidSc9+0NIn/2uoea1yOScQ0C\nahJptt65AnUyX69FiEQrlL9IN/bb/fu79Hkkv1/w1pzWRyATSL/RbD2ShlJX5xqZRDLfGe+z\nRZr25pae2eCXnNacNqxP5A+/0fJNUjrIuUYWkXabFWk/jUgLX9jnlZuanNZkumutVGoUyb0G\ndTI3477dNkX6Yvv2SqSf/Yq7dokiuVdLpIbrZKeari7SZZsiNafheN3S1yPFiOQ3IKnrJVLB\nvp3SQe5V8oi026xIzU+r0ulbuyynSD5HG9SM1uS5x4qJ1CeSxyrkIjX93YSMz7P4dJFsrCyS\nPaGta5NJpNlutk6JR66oxyoQyYMtieS3ahKU+3ZUcbRBCxBpt1mRvrtdO+29T5YSyfS9axJa\nzXLPdZNIH5IMgSMiWaKHHGvIIFJzGe4UWYRIJEIFiHQcDjYsPbMht0ikJqWJZA0dGMu6BS+P\nciVzM4ikuQnKFkS6t4e/GXsd2V2ztGaRKE2K37fzCB7Ub9atrC3Sn+FuQlsQ6dCekGWsebGD\nZum6ImnSOSgNfTXxIHZI8ose1nXmrZQiknzvhi2I1M0Mav+3+BShDYjkHT609/TbmcZNa+Hs\nIqk3QdmSSK/Fpwi5D9s5stmdg96J7GQXY1LIBoI7ULMhv0MNOUUahiT5kvMtiHRufxut8xvJ\nOSRpsjkwA0NS2U64SIEbiOhDeUuliKTcu2ELIj33rL+6b/EpQjQiOTYTmM1mQkWK2EREL4rb\n8tqzyy2S7iYoq4o0fRBDwOHv56kVaY3zSC6RXMnslX0RCa1F3naOzUZ1ZFOUSPqboGxCJCtF\niRS+Z6eLEon8IynPb7OorlSqaC+eW6QLRFL5AJEITfLdcsJWYrpSqaG9OEQKYO2nmrshEsln\nU1HprHymbp1scz5RPZrnNyBlFonbt4NI3NLI7btZUKS442earXuZFLgx0+b9USroKJ9dJPUm\nKFsQaeBxPL00Hy8mks9lrsuIZFxLs/nETTmrENQ43wmrmUWa78u1SZGWnyLkvFpOya9okWKn\nGARv3nszQdXwbFo5Iin3btiSSMs/aCxUpASPYue8Obef7T5g4Q1z7tmZMo9apOl2QoWJFCVU\nhEg33XUUHyOSV5bb19L8SMp5Y8rQdjkHpAVFumxZJC0riqRkVppIHonuWkt/oXs2AltVjkjq\nJecQ6YNEcprkXGlhkcIm9Lkvjl1YpB1E4ihcpLBNRmWtWIGFTQqYAOU+1LCkSNIl5xBpOZHc\nB8CSRbKa5LPO8iL5T+krSSTlknOIVJBIqXt2+qjuSI4qZMevOR63PVlUJPFKWYiU9bEutn07\nNaFIRDKp5LfCCvt2vnWba2YsvKBI8pWyEKlskeK2HJar6w9JXrudHmdjlxVJuFIWIq0lkiab\nyESSooeUXkkk9xFFjz27RUWSrpSFSJ8pUhjafF0cW7X8bp6/sEj9uaRLeSIFCVW5SLpMWk2k\nIoYkubXSMotIzsTKIpJ4pSxEWlAkPgU0abSeR4UMST36/rHt2a0m0gUizSwpkn3e2ooiFWWS\nFtue3aoiDdclQaRFRZqywJ4rSvHsOOqxPgWKNF3gd4FIzdIijWlgzxW59AKI9UgxaSeQEEgO\nOrxU676WSNMFfhf60BBp5t/VaJIpWXRll8FekZB81xEZTlclTd1XFOkKkUbyXmrOZ9GYVc2Q\ntlKWKVnHle/I+56rz1AJU330700OzYXC4onvx2pp6z/kzby8/4ArL7z/U983wvs/9X0jvP8T\n3neXUzTC8j/1fSO8/1PfN8L7P0v7nO1V2m/7/msWSZs4u52xfE/e93x9doGJvxOQl48f+8dT\n3gvrK/UvQKQrRFp61+6N757dkrt27itlTew47EW8Q+oCjC81VXfu6uTatZuvS1IfmES97aTQ\nNioWyZks5pzJiL0ullr7WJo0hpsAAA1sSURBVJKo0ryqruarinSBSD0QaUDOe886hxX17AvN\nusNLXc1XFGm6Lkl5YBJEIiNMpJX37JoIkwKHmXiVDHt23omVUaTxcgrlgUnU24ZI4cmiS5oF\ncFVHrXCoGLEmlS3ScGGSdDNw6m1DpPBk0eTMEig576huhBVxg5Jpz847sXKKNM8UgkiZ+Dck\ngXe2SJ9kq5gJU/bqUzt6cAldj19DqK93YuUWabieQrj1alEiTR9ogEjkeGd86oGD4B3C6bVQ\nX+/EyipSa9J1NAkiZSBIpPV/IjWaO4Fr605wUihkfeOAVJhInUnXbNuGSLEiZauXGUsKx2pg\naGtADOOAVIpI7Szw0ST1fFJRIk0LOCASPbqM1ymQPAk1IIx5QCpGpHby6mCSej6paJHaToVI\nGbBlse8kBj98QwmFpNr6502+lOz/vY531VfPJ5Uq0tSpnyVSGR5p71kp49MczyZ7xBMKSLX1\nz5t8KTn8e73OJ2Z1D5dN3zatSFwPQ6QcGNOd2CIxtr2IsUv88yZfSo7/DiJ1J2Y1TyBzhG6H\nMscpXcJaX/tKXvperUQkT5OKFemayyIxuKWAuUu8E2sZkaYTsxf1eS+uNQeRLv3OYc4Hpo+b\ngEhZCfOAAqtK1gGpKJGai2CSfFNw/ZqcQFdBJI2KRLW+bE2kbLVyEG5COkaVdtYBqSyR/gST\n5nHFtOaczAaRpAgktZY30fJRIpUyIK0jkl4l5aOpjsGJtYhI/R/7waQr74myplzAKBLx7SCU\nTbRApDxEqpCK9EtMs8M3V7FQkXo/umr7e+JRgKrWum1fIVIu0nxIQD7MLg9RcxWLFukyVJ1O\nJJprBvWhr58uUrZKOUn1IQHroXauiqWKNPzRnxrg54nyx0NnmnelAh291nLzk/pEWtOkDsOR\nB66G5Yo0nKKZVbpaklk3BPNrijp4TCwPHez6LX2SSAXt2TXrm6SHq2DBIk3JzCmhS2arRMY1\nLZUy7BQaRZq3UotIPiaVJVKRJvH1K1mk4TC4mKtW5GxXRibHmdpGHYlcIvFbr+Y3UpRI2erk\nBWXSU8STY5YtEn/O03/AuRqHLOuZWtuemyE0H7kFIuWDKOFpIqphCxdJ+sGvl8hn/0tWLmDP\nzViAj9l3LUTKSHqyE8XUxi5dpD9hXmhEtqvDB41IopmfJ1JxHoUlfZ6olk1UINL8Wyn1PJJm\nDIsTSTsaXiFSXmJznC6uZSN1iOSYkyqdJmosx97kcSRYJOmXGURakJgEJwxs3UwtIrX/WrPd\nGFrWwaCBFFp/TEOy6INF2iklslXJn/D0pgxt31BNIkWHlnyxHaZwC2TcObxWJJLTpBI9cqd7\nztiGTUXnbI0iyb+y1HO8nv7sNLOOINKCZNPIHd2wra2JpJ4mCvTF/SPqk0Qqcs+usaZ69g3o\nN7Y9kf76GxTx6R804GxNpBI9akyZnn0DEEkqoGgw+JN0ZL0P3VRzg8iaRdJlev4tmDe4VZH+\nwuek2gvMoZvPEanUPbuBBWrm69GWRfrjD+MliSSGbioTyX7ntoI9WgaI5Bc6XSQ5dFOTSNYh\nqfABaTE8PIJIwRfAigU0oZtPEsmWO1vC5RFE6v/VeiKv6SzwV6dIfrcS3bRHjaKSvDg6JT9L\nJMrQTVUiWYYk7NkpmPshOm8gkqlA80EiOf4Mg/S8gUimAk11IhlMgkhu0vMGIpkKNHWJZByS\nsGfnQXreQCRTgeZzRMKA5CI9byCSqUBTn0g6kzQfZ6tNvaTnDUQyFWgqE8kwJMEjH9LzBiKZ\nCjRViNTYlIFIvqTnDUQyFWgqFEk1CXt2XqTnDUQyFWhqE0k7JMEjL9LzBiKZCjQQ6fOhyxuI\nZCrQ1CGSfd8Oe3Z26PIGIpkKNNWJpBl+4JEduryBSKYCzQeIpNnZy1aVKqHLG4hkKtBUIpJl\n3053HC9bVaqELm8gkqlAU59I0gj08R7ZvrvINaPzBiKZCjS1i6Q9r5StJkb4DqUM2Xh8uc5K\n0eUNRDIVaOoUace/KcCjBJHovlwxJGlowozzK1Bd6KYWkQxDknbGULaKqJg6ODBCLdkOkUwF\nmkpFGgRafcfO1MHxa6Z/uXWG/iNac6XQTY0iCY8PWHVAgkiZClQXuqlGJINJ63gU1MGRa0YX\nqDP0H9GaK4Vu6hRJ2L+DSORpsULoP6I1VwrdVCtS55LGo8JEyp0WmgJ1hv4jWnOl0E3NImnJ\nVgueIr47U4E6Q/8RrblS6KYekfI91tifor47U4E6Q/8RrblS6ObTRMpWiY6ivjtTgTpD/xGt\nuVLoBiKFUNR3ZypQZ+g/ojVXCt1UJJKPSZlqUOR3ZyrgXvO/ljyhiQpUF7qBSD4U+d2ZCkCk\nFUI3nyUS9ZaL/u5MBYwL/hsFEv79DyJBpEwi1fHdmQpApBVCNzWJ5DSJaoN1fHemAhBphdDN\nJ4mUvKG6vjtTgXCRZKEgUnCB5oNEit9And+dqYCyQPYFItGHbj5HpPC4lX93pgIQaYXQTVUi\n2UwKiPch352pQLxIUwHv0BX2GERqIRmPPuS7MxUgEMlUACKZCjSViWQyySvOh313pgLTgmBP\n3CIJQtXZYxCpJ2E8+rDvzlRgCZHGDyrsMYg0EO7Rh353pgIQaYXQSprGs5RIikmm1YroYOOC\nfCkZ70m4SHLty872jKHVNI2GUiTG+GhyDT0HoyI62LggX0quIdJUIGOPFR1aTFMxfwMhFImJ\n4QJVL6qD5Q/Sk5kg2xcITddjGb8MwtBCmkr5GwidSEyK19dQV3eBNTtYTihj6KKyfYHQa3wZ\na4RuOJHk/A0kt0iJX42zgM+akYlFmswZsn2N0N5fhjF0xu85vEBTukj/AKgHJX+j05+CxIpw\nUB1IqS1kFZWsI2RERIj0MSGrqGQdISFSSxVfVRnfPkKSRYRIHxOyikrWERIiAbAOxYiUeEIL\ngHUp5YRs4hQLAFamlClCAGwXiAQAARAJAALWF2naM2XzK+ljopBTa+lCTu/CQ2oaSNluLiRV\nu8VaxrabX4Wo4YaIKe0OZnWRpmMl80ETpv84NWTboeLyNWu51Xbzqwh/kuJDGiKmtDuctUWa\nmsp98Uz/cWLI9n9MWL5mLcVV1UDJlRRqStNu2lpO9SLrSjFiSrsjKEYk7i3dV8WH1KVYasj4\nWoprUKeo8I6s3WovFC1SQrsjKFAkJnzMmtBKGkKKi8hCRtZSiChXi0AkjQK0IdNrybQNT+lK\npUax33cEZYnEjx7T29BfisYcmlOCMOS4LDCkpoGEIvEhpRpThZS/pvCQBpESutIg0iYONjhF\nSh0+NCKRhpTrnBZRajdFSHERWUiuznEh1eamjkiq2lsdkZj4v/S/zMYvf/WQ6p93OVD7ZzQh\n6+fOmz8sKaT0XSd84YaIauislCQSn0Xjy8SsVzOBOKR2M4ERG61ISZXkQkpxiEIypVRgSDZA\nKJIYcVq0OZHUP/jJIvEh1dgEIZmyPDSi+IKk3dKLhqDdcsjodourkI5IUkQ1dFbWFqmrgCFD\nmfAiOWTD9ylNSKYsD47IhheE7ZZCNgTt1tQyNaT0OqHhhohNQrvDWV2k8ZDKOC4z8eMm5pCL\nIeT8pRGFnPYlYkJKDZTbzeIrKYdsGuKQc7uTQkoKkE0Rovi+g1lfJAA+AIgEAAEQCQACIBIA\nBEAkAAiASAAQAJEAIAAiAUAARAKAAIgEAAEQCQACIBIABEAkAAiASAAQAJEAIAAiAUAARAKA\nAIgEAAEQCQACIBIABEAkAAiASAAQAJEAIAAiAUAARAKAAIgEAAEQCQACIBIABEAkAAiASAAQ\nAJGWoX8QyuH2khf8OFd9nRm7jVFMS7RhF3iaCRhAVy/D+FzG/VP8/OD+Ak7v1b7GKKYlEgem\nKQ4ygq5ehj6nn0d21H3uWPUZsSSgboAAdPgyjJl9YD/azz1WpVkC8oAOX4Yxs3/Yuf3/e6ds\nf2uGPb728/uB7e/8Cu8PDvexBOOivIehE9t/cUu4dW97dnxOYafQXaR5TUAPRFqG0YUXOzTN\nV+/AbRbp1L3gdvuO4wcakfas+200LuHW7dbavwSRpkjzmoAeiLQMkgvfTfPNpseFv8ep46t5\nHefdvm+2/21+9205fjetX/ld9t762C/h1v1uX557QcfiXCRuTUANRFoGQSTuVf/2xNrD4i92\nGpedOqd+hoFEiMLYQ1iZW/fULnqxPS8SF4lbE1CDXl0GSaTnz9eRE2k8OM6k4lyR6f3syPi/\naV15fe4zaU1ADXp1Gcb0fXaDzHHSBiJ9COjVZRjT97v9BXNmh/vPUxBJX9xPJGUjEGl50KvL\nMJ9HegxveJFO0tml6ZfNqXGKxK17tPxGOkGkrKBXl0GY2dD+6v8dfyO1kxO6Q2vNfT7YYD1q\nN8Xr/sete28Py936o3bPRnPUrmkgUibQq8sgzLW7DW8e7QjVDiDDjyZuIt589sclEr/ueB5p\nCKucR5LCATrQq8vQm3McToae3y8f3e7W49CJ1E4/YGd+5tx9P8xHcIrEr/tW9NS+6sMOMxv2\n08wGKRygA70KAAEQCQACIBIABEAkAAiASAAQAJEAIAAiAUAARAKAAIgEAAEQCQACIBIABEAk\nAAj4H/J6ZPcKDkFlAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 300,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "scaling_parameter=max(result_withcut$total)/max(result_withcut$upper[!is.na(result_withcut$upper)])\n",
    "range = c(as.Date(\"2020-02-01\"), (datestar-17))\n",
    "options(warn=-1)\n",
    "\n",
    "figure_data %>% \n",
    "    ggplot() + \n",
    "    geom_bar(aes(x=onset, y=value, fill=subject, group=subject),stat='identity', width=0.7) +\n",
    "    scale_fill_manual(values=c(\"#FAAB18\",\"gray35\")) +\n",
    "    geom_line(data=result[!is.na(result$Rt),],aes(x=onset,y=Rt*scaling_parameter),color=\"#1380A1\",size=1) +\n",
    "    geom_ribbon(data=result,aes(ymax=result$upper*scaling_parameter, ymin=result$lower*scaling_parameter, x=onset), \n",
    "                fill=\"#1380A1\", alpha = 0.4) +\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",
    "    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": 16,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "write.csv(result, '../results/rtmle.csv')\n",
    "write.csv(result_withcut, '../results/rtwcmle.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}