kenmomd
diff --git a/chiba_kanagawa_saitama_onset_dates.csv b/chiba_kanagawa_saitama_onset_dates.csv
diff --git a/chiba_kanagawa_saitama_rt_by_onset_master.csv b/chiba_kanagawa_saitama_rt_by_onset_master.csv
diff --git a/covid-19_chiba_kanagawa_saitama_rt.R b/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)
	dates	counts
	2020-02-16	1
	2020-02-17	0
	2020-02-18	0
	2020-02-19	3
	2020-02-20	0
	2020-02-21	1
	2020-02-22	1
	2020-02-23	1
	2020-02-24	4
	2020-02-25	2
	2020-02-26	3
	2020-02-27	4
	2020-02-28	7
	2020-02-29	5
	2020-03-01	6
	2020-03-02	6
	2020-03-03	5
	2020-03-04	8
	2020-03-05	4
	2020-03-06	1
	2020-03-07	4
	2020-03-08	5
	2020-03-09	9
	2020-03-10	3
	2020-03-11	5
	2020-03-12	5
	2020-03-13	11
	2020-03-14	5
	2020-03-15	5
	2020-03-16	10
	2020-03-17	18
	2020-03-18	17
	2020-03-19	14
	2020-03-20	19
	2020-03-21	16
	2020-03-22	23
	2020-03-23	35
	2020-03-24	34
	2020-03-25	43
	2020-03-26	36
	2020-03-27	79
	2020-03-28	64
	2020-03-29	66
	2020-03-30	102
	2020-03-31	88
	2020-04-01	105
	2020-04-02	88
	2020-04-03	99
	2020-04-04	91
	2020-04-05	92
	2020-04-06	117
	2020-04-07	93
	2020-04-08	90
	2020-04-09	65
	2020-04-10	68
	2020-04-11	60
	2020-04-12	57
	2020-04-13	62
	2020-04-14	68
	2020-04-15	51
	2020-04-16	38
	2020-04-17	40
	2020-04-18	32
	2020-04-19	35
	2020-04-20	29
	2020-04-21	29
	2020-04-22	36
	2020-04-23	20
	2020-04-24	34
	2020-04-25	27
	2020-04-26	24
	2020-04-27	28
	2020-04-28	25
	2020-04-29	15
	2020-04-30	22
	2020-05-01	21
	2020-05-02	21
	2020-05-03	8
	2020-05-04	13
	2020-05-05	10
	2020-05-06	9
	2020-05-07	14
	2020-05-08	10
	2020-05-09	7
	2020-05-10	11
	2020-05-11	8
	2020-05-12	11
	2020-05-13	8
	2020-05-14	8
	2020-05-15	10
	2020-05-16	6
	2020-05-17	9
	2020-05-18	2
	2020-05-19	5
	2020-05-20	3
	2020-05-21	0
	2020-05-22	3
	2020-05-23	3
	2020-05-24	5
	2020-05-25	2
	2020-05-26	5
	2020-05-27	4
	2020-05-28	4
	2020-05-29	4
	2020-05-30	1
	2020-05-31	4
	2020-06-01	2
	2020-06-02	6
	2020-06-03	0
	2020-06-04	2
	2020-06-05	7
	2020-06-06	4
	2020-06-07	5
	2020-06-08	5
	2020-06-09	6
	2020-06-10	2
	2020-06-11	7
	2020-06-12	2
	2020-06-13	3
	2020-06-14	7
	2020-06-15	1
	2020-06-16	7
	2020-06-17	5
	2020-06-18	9
	2020-06-19	7
	2020-06-20	11
	2020-06-21	7
	2020-06-22	16
	2020-06-23	16
	2020-06-24	16
	2020-06-25	13
	2020-06-26	12
	2020-06-27	19
	2020-06-28	14
	2020-06-29	18
	2020-06-30	17
	2020-07-01	18
	2020-07-02	29
	2020-07-03	23
	2020-07-04	25
	2020-07-05	32
	2020-07-06	40
	2020-07-07	25
	2020-07-08	40
	2020-07-09	40
	2020-07-10	36
	2020-07-11	42
	2020-07-12	39
	2020-07-13	62
	2020-07-14	43
	2020-07-15	50
	2020-07-16	42
	2020-07-17	53
	2020-07-18	39
	2020-07-19	55
	2020-07-20	56
	2020-07-21	65
	2020-07-22	70
	2020-07-23	64
	2020-07-24	80
	2020-07-25	79
	2020-07-26	75
	2020-07-27	109
	2020-07-28	94
	2020-07-29	114
	2020-07-30	99
	2020-07-31	84
	2020-08-01	105
	2020-08-02	75
	2020-08-03	117
	2020-08-04	101
	2020-08-05	95
	2020-08-06	93
	2020-08-07	69
	2020-08-08	94
	2020-08-09	80
	2020-08-10	92
	2020-08-11	96
	2020-08-12	102
	2020-08-13	118
	2020-08-14	96
	2020-08-15	78
	2020-08-16	75
	2020-08-17	93
	2020-08-18	85
	2020-08-19	83
	2020-08-20	56
	2020-08-21	72
	2020-08-22	63
	2020-08-23	65
	2020-08-24	74
	2020-08-25	72
	2020-08-26	76
	2020-08-27	76
	2020-08-28	77
	2020-08-29	47
	2020-08-30	67
	2020-08-31	73
	2020-09-01	74
	2020-09-02	63
	2020-09-03	65
	2020-09-04	59
	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)