chiba_kanagawa_saitama_onset_dates.csv

chiba_kanagawa_saitama_rt_by_onset_master.csv

covid-19_chiba_kanagawa_saitama_rt.R

if (!requireNamespace("pacman", quietly = T)) install.packages("pacman")
pacman::p_load(tidyverse, lubridate, scales, glue, patchwork, ggthemes, incidence, EpiEstim, curl)
# 「含まれるか」を意味する%in%のnot equal版
'%!in%' <- function(x,y)!('%in%'(x,y))

#######################################################
#jagjapan　首都3県の実効再生産数---------------------------
#######################################################
URL = "https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv"
PATH <- "jagjapan.csv"
download.file(URL,PATH)

jagjapan_master <- read_csv(PATH) %>% dplyr::select(6:16)
#日付
jagjapan_master$確定日 <- lubridate::mdy(jagjapan_master$確定日)
jagjapan_master$発症日 <- lubridate::mdy(jagjapan_master$発症日)
#TimeLag
jagjapan_master$timelag <- jag$確定日 - jag$発症日

#受診都道府県が"千葉県","埼玉県","神奈川県"かつ確定日が3月以降の陽性者を抽出
pref <- c("千葉県","埼玉県","神奈川県")
confirmed_0301 <- jagjapan_master %>% dplyr::filter(確定日>="2020-03-01" & 受診都道府県 %in% pref)
#居住都道府県で抽出
#residence_0301 <- jagjapan_master %>% dplyr::filter(確定日>="2020-03-01" & 居住都道府県 %in% pref)

#発症日がNAの例数
confirmed_0301 %>% dplyr::filter(is.na(発症日)) %>% count()
#n=5398

#NA削除
confirmed_0301 <- confirmed_0301 %>% dplyr::filter(!is.na(発症日))

#incidenceオブジェクト作成
dat_date <- confirmed_0301$発症日
#確認
new_i <- incidence(dat_date)
new_i
plot(new_i)

date_onset_latest <- length(new_i$dates)

#データが不十分な直近をカット
plot(new_i[12:date_onset_latest-12])
T <- nrow(new_i[12:date_onset_latest-12]$counts)
#7日間windowでRtを算出する
t_start <- seq(2, T-6)
t_end <- t_start + 6

#Rt算出
estimated_Rt <- estimate_R(new_i[12:date_onset_latest-12],
                           method = "uncertain_si",
                           config = make_config(list(
                             mean_si = 4.8, std_mean_si = 0.7109351,
                             min_mean_si = 3.2, max_mean_si = 6,
                             std_si = 2.3, std_std_si = 0.7695178,
                             min_std_si = 0.5, max_std_si = 2.5,
                             t_start = t_start,
                             t_end = t_end)))

#Rtをcsvに保存
chiba_kanagawa_saitama_rt_by_onset <- as.data.frame(estimated_Rt$R)
write.csv(chiba_kanagawa_saitama_rt_by_onset, file = "chiba_kanagawa_saitama_rt_by_onset.csv", row.names = FALSE)

chiba_kanagawa_saitama_onset_dates <- cbind.data.frame(estimated_Rt$dates, estimated_Rt$I) %>% 
  as_tibble() %>% 
  rename(dates=`estimated_Rt$dates`,counts=`estimated_Rt$I`)
#発症日をcsvに保存
write.csv(chiba_kanagawa_saitama_onset_dates, file = "chiba_kanagawa_saitama_onset_dates.csv", row.names = FALSE)

#全体のグラフの確認
plot(estimated_Rt, legend = FALSE)

#描画範囲
lims <- as.Date(strptime(c("2020-03-01","2020-09-05"), format = "%Y-%m-%d"))
#今日の日付
today <- lubridate::today()

#新規陽性者数グラフ
p1 <- plot(estimated_Rt, "incid") +
  theme_clean()+
  labs(x="", title="COVI-19 千葉県・神奈川県・埼玉県の発症日別陽性者数" ,y="")+
  scale_x_date(limits = lims, date_breaks  = "3 days", labels = date_format(format = "%B%d日"))+
  theme(axis.text.x = element_text(angle = 90,hjust = 1)) +
  theme(legend.position = "none")+
  scale_y_continuous(expand = c(0,0), breaks = seq(0,500, by=10))+
  scale_fill_manual(values = c("#AD002A"))
p1

#Rtグラフ          
p2 <- plot(estimated_Rt, "R")+
  theme_clean()+
  theme(legend.position = "none")+
  labs(title = "発症日ベースの実効再生産数", x="", y=bquote("Estimated" ~ R[t]), caption = glue("mean serial interval = 4.8 days, SD = 2.3 days\n（ヽ´ん`）@kenmomd {today}作成 \nDATA：J.A.G JAPAN"))+
  geom_hline(yintercept = 1, colour = "steelblue") +
  scale_x_date(limits = lims, date_breaks  = "3 days", labels = date_format(format = "%B%d日"))+
  theme(axis.text.x = element_text(angle = 90,hjust = 1)) +
  scale_y_continuous(expand = c(0,0), breaks = seq(0,6, by=0.2))+
  coord_cartesian(ylim = c(0.2,2.5)) #y軸範囲
p2

pp <- p1/p2
pp

#グラフ画像の保存
ggsave(filename = glue("Chiba_Kanagawa_Saitama_COVID-19_Rt_by_onset_{today}.png"), plot=pp, width=15, height=10, dpi = 160)

#2軸プロット
chiba_kanagawa_saitama_rt_by_onset <- read_csv("chiba_kanagawa_saitama_rt_by_onset.csv")

p3 <- chiba_kanagawa_saitama_rt_by_onset %>% 
  select(c(dates,incidence,"Median(R)","Quantile.0.025(R)" ,"Quantile.0.975(R)")) %>%
  mutate(rescaled_Rt = (`Median(R)`)*60-60) %>%
  mutate(rescaled_Rmin = (`Quantile.0.025(R)`)*60-60) %>%
  mutate(rescaled_Rmax = (`Quantile.0.975(R)`)*60-60) %>%
  pivot_longer(cols = c(-dates,-rescaled_Rt,-rescaled_Rmin,-rescaled_Rmax), names_to="var", values_to="value") %>%
  #dplyr::filter(var != "incidence") %>% 
  dplyr::filter(var %!in% c("Quantile.0.025(R)", "Quantile.0.975(R)","Median(R)")) %>% 
  ggplot(aes(x=dates))+
  theme_minimal(base_size = 13, base_family = "HiraginoSans-W4")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(axis.text.y  = element_text(angle = 0))+
  scale_x_date(date_breaks  = "1 week", labels = date_format(format = "%Y-%m-%d"))+
  scale_y_continuous(sec.axis = sec_axis(~./60+1,
                                         name = bquote("Estimated" ~ R[t]), breaks = seq(0,7,0.2)),  breaks = seq(-100,300,20))+
  coord_cartesian(ylim = c(-40,120))+
  geom_bar(aes(y = value), stat="identity", alpha=0.7, fill="orange",colour="orange", width=0.2)+
  geom_hline(yintercept=0)+
  geom_ribbon(aes(ymin = rescaled_Rmin, ymax = rescaled_Rmax), alpha = 0.4, fill="light blue", colour="light blue") +
  geom_line(aes(y=rescaled_Rt), size=0.8, colour="steel blue")+
  labs(title="千葉県・神奈川県・埼玉県の発症日別COVID-19陽性者数と発症日基準の実効再生産数",
       x = "", y = "Cases", 
       caption = glue("mean serial interval = 4.8日, SD = 2.3日, 7日間window\n（ヽ´ん`）@kenmomd {today}作成 \nDATA：J.A.G JAPAN")) 
p3
ggsave(filename = glue("chiba_kanagawa_saitama_Rt_by_onset2_{today}.png"), plot=p2, width=13, height=10, dpi = 160)