updated

DATA/CORONA_TR.csv

@@ -84,4 +84,7 @@
 01.06.20,164769,827,4563,23,651,974,31525
 02.06.20,165555,786,4585,22,633,974,32325
 03.06.20,166422,867,4609,24,612,931,52305
-04.06.20,167410,988,4630,21,602,926,54234
+04.06.20,167410,988,4630,21,602,926,54234
+05.06.20,168340,930,4648,18,591,1622,36846
+06.06.20,169218,878,4669,21,591,1922,35846
+07.06.20,170132,914,4692,23,613,2647,35335

OUT_07_Jun_2020/CSVS/parameters_fit_mrelax_0.csv

@@ -0,0 +1,11 @@
+"mean" "se_mean" "sd" "2.5%" "25%" "50%" "75%" "97.5%" "n_eff" "Rhat"
+"R0" 4.14913022198168 0.0171412700679499 0.536464526145957 3.05476850387226 3.77227112315959 4.2122332433367 4.58356997347662 4.94629017514915 979.481017938959 0.999134708925276
+"r_lock_1" 0.202230507243658 0.00118714148238749 0.0357842267048068 0.1557015845329 0.173894624241584 0.194293800819742 0.223364936534366 0.28527649838917 908.611672025837 0.999948700298437
+"gamma_s" 0.376182261883525 0.00327098836815715 0.118987961410773 0.202098440599599 0.288819185959371 0.354724871702634 0.441144827014135 0.663482607401664 1323.26872547164 1.00037538018373
+"gamma_H" 0.164475263497818 0.00111258805833209 0.0286655068301453 0.12283212546864 0.14434575439295 0.15892713168163 0.179225891349763 0.233233777376614 663.820196060723 1.00713290444574
+"gamma_ICU" 0.163073834540889 0.000914950734895629 0.0353682409512046 0.109616175582454 0.138672894876803 0.158175267458254 0.180405885414893 0.246849695985906 1494.27833763923 1.000725004756
+"eps_H_ICU" 0.199314011794776 0.00222542564682253 0.0660567647144907 0.107762682143303 0.153259869735022 0.18765662965471 0.231603115935798 0.363604109404325 881.065956904847 1.00014703046975
+"eps_H_x" 0.0483125243484855 0.000434164509779369 0.0138201629772885 0.0268884818915933 0.0385815590491009 0.0466514730231838 0.0559359698872961 0.0797727065182705 1013.25251344432 1.00054759484558
+"eps_ICU_x" 0.0996596228056717 0.00160630437722676 0.051058149682627 0.0185099235486842 0.0603605455558992 0.0934979206038248 0.133182615224749 0.211167993799205 1010.35607871321 1.00437778733785
+"r_d_s" 0.0680956311731065 0.00148722963557997 0.0449334367719451 0.0135915959209827 0.0358392735991378 0.0582429341581251 0.0878231203067855 0.181506834736344 912.815934442016 1.0001684220948
+"r_d_a" 0.178845164696777 0.00275985886592931 0.102288677445906 0.0337770279303058 0.100197450671172 0.161273216428045 0.236038489599487 0.426583808279533 1373.66672904448 1.00045661359858

analysis_plots_RETRO.R

@@ -0,0 +1,112 @@
+# Setup
+rm(list=ls())
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+path2save = paste0("OUT_",format(Sys.time(), "%d_%b_%Y"));
+dir.create(path2save)
+dir.create(paste0(path2save,"/CSVS/"))
+dir.create(paste0(path2save,"/FIGS/"))
+dir.create(paste0(path2save,"/RDATA/"))
+
+args <- commandArgs(TRUE)
+print(args)
+
+m_relax_in = 0 # slope of relaxation -> if 0, no relaxation
+
+library(rstan)
+library(zoo)
+library(Rcpp)
+options(mc.cores = parallel::detectCores())
+rstan_options(auto_write = TRUE)
+source("prepare_data.R")
+source("setup.R")
+
+days2add     = 45; #ADDITIONAL DAYS FOR SIMULATION
+date_simul   = date_end + days2add;
+
+data_list = list(
+  pop_t=pop_size,
+  tswitch_1=as.numeric(date_control_3-date_data+1),
+  trelax=as.numeric(date_relax-date_data+1),
+  r_end=1,
+  m_relax=m_relax_in/100, 
+  K   = test_fit_vec_1[1],
+  mu  = test_fit_vec_1[2],
+  sig = test_fit_vec_1[3],
+  D=as.numeric(date_end-date_data+1),
+  k_daily_cases  = daily_cases_data,
+  k_icu          = icu_data,
+  k_daily_deaths = daily_deaths_data,
+
+  p_pi        = c(1,999),
+  p_R0        = c(3,1),
+  p_tau       = c(2,3.2),
+  p_gamma_s   = c(2,3.2),
+  p_gamma_H   = c(2,3),
+  p_gamma_ICU = c(2,3),
+  p_eps_H_ICU = c(2,3),
+  p_eps_H_x   = c(2,10),
+  p_eps_ICU_x = c(2,10),
+  p_r_d_s     = c(2,10),
+  p_r_d_a     = c(2,10),
+  p_r_lock_1  = c(1,1),
+  p_phi       = 1/100,
+
+  t0=0,
+  t_data=1,
+  S=as.numeric(date_simul-date_data+1),
+  E=days2add, 
+  ts=1:as.numeric(date_end-date_data+1),
+  ts_pred=1:as.numeric(date_simul-date_data+1),
+
+  r_c=0 #dummy control
+)
+
+load("/Users/burcu/Dropbox/corona-tr-modeling/OUT_31_May_2020/RDATA/T_modelTR_mrelax_0.RData")
+T_modelTR_old = T_modelTR
+model_end = ymd("2020-05-31")
+D_old=as.numeric(model_end-date_data+1)+days2add
+
+pp = c("predicted_daily_cases","predicted_daily_deaths","predicted_current_icu")
+model_output_fitted_old = summary(T_modelTR_old,pp)[[1]] %>%
+  tbl_df() %>%
+  mutate(t=rep(1:D_old,3),
+         date=date_data+t-1,
+         eta="-100%",
+         type=rep(pp,each=D_old)) %>%
+  mutate(populasyon=factor(type,levels=pp,
+                           labels=c("Gunluk Vaka Sayisi","Gunluk Vefat Sayisi","Yogun Bakim Hasta Sayisi")))
+
+D_data = as.numeric(date_end-date_data+1)
+populasyon_data=factor(rep(pp,each=D_data),levels=pp,
+                  labels=c("Gunluk Vaka Sayisi","Gunluk Vefat Sayisi","Yogun Bakim Hasta Sayisi"))
+
+agg_data         = tibble(allDates_agg,daily_cases_data,daily_deaths_data,icu_data);
+observed_data    = c(daily_cases_data,daily_deaths_data,icu_data);
+observed_dates   = c(allDates_agg,allDates_agg,allDates_agg);
+agg_data_all     = tibble(all_dates=observed_dates,all_data=observed_data,populasyon=populasyon_data);
+
+
+old_data_length     = as.numeric(model_end-date_data+1);
+populasyon_data_old =factor(rep(pp,each=old_data_length),levels=pp,
+                       labels=c("Gunluk Vaka Sayisi","Gunluk Vefat Sayisi","Yogun Bakim Hasta Sayisi"))
+observed_data_old   = c(daily_cases_data[1:old_data_length],daily_deaths_data[1:old_data_length],icu_data[1:old_data_length]);
+observed_dates_old  = c(allDates_agg[1:old_data_length],allDates_agg[1:old_data_length],allDates_agg[1:old_data_length]);
+agg_data_all_old    = tibble(all_dates=observed_dates_old,all_data=observed_data_old,populasyon=populasyon_data_old);
+
+# # 
+ggplot() +
+  geom_ribbon(data=model_output_fitted_old,aes(x=date,ymin=`2.5%`,ymax=`97.5%`,fill=populasyon),alpha=.5) +
+  geom_line(data=model_output_fitted_old,aes(x=date,y=`50%`),colour="black") +
+  facet_wrap(~ populasyon ,scales="free",nrow=2) +
+  geom_vline(xintercept=date_control_3,linetype=2) +
+  geom_vline(xintercept=model_end,linetype=1) +
+  geom_vline(xintercept=date_end,linetype=1) +
+  scale_colour_manual(values=c("grey20","black"),guide=FALSE) +
+  scale_alpha_manual(values=c(1,0),guide=FALSE) +
+  scale_x_date(date_breaks="2 weeks",date_labels="%b %d") +
+  theme(axis.text.x=element_text(angle=45,hjust=1)) +
+  geom_point(data=agg_data_all,aes(x=all_dates,y=all_data,fill=populasyon),shape=21,size=1,colour = "black", fill = "black") +
+  geom_point(data=agg_data_all_old,aes(x=all_dates,y=all_data,fill=populasyon),shape=21,size=1,colour = "black", fill = "white") +
+  labs(x="Gun",y="Birey Sayisi")
+# scale_y_continuous(trans = 'log10')
+ggsave(paste0(path2save,"/FIGS/figure_fit_mrelax_",m_relax_in,"_RETRO.png"),width = 10, height = 5)

model_fitting.R

@@ -1,6 +1,6 @@
 # Setup
 rm(list=ls())
-setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+setwd(getwd())
 path2save = paste0("OUT_",format(Sys.time(), "%d_%b_%Y"));
 dir.create(path2save)
 dir.create(paste0(path2save,"/CSVS/"))
@@ -61,9 +61,9 @@ data_list = list(
   r_c=0 #dummy control
 )
 # # IF .rds NOT compiled (run in case of change in model)
-M_model_TR     = stan_model("MODELS/model_TR.stan")
+# M_model_TR     = stan_model("MODELS/model_TR.stan")
 # # IF .rds  compiled 
-# M_model_TR = readRDS("MODELS/model_TR.rds")
+M_model_TR = readRDS("MODELS/model_TR.rds")
 ####### FITTING - DEBUG MODE
 # T_modelTR      = sampling(M_model_TR,data = data_list,iter=5,chains=1,init="random") 
 ####### FITTING - SHORT VERSION

prepare_data.R

@@ -1,5 +1,5 @@
 # Setup
-setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+setwd(getwd())
 source("setup.R")
 library(anytime)
 library(readr)