# Plotting Rt for covid-19 in Tokyo in severl methods. # July 3, 2021: Made from libusejapant.R # Instantaneous Rt estimtion series for Tokyo if (!require(COVID19)) { install.packages("COVID19", dep=TRUE) library(COVID19) } if(!require(EpiEstim)) { install.packages("EpiEstim", dep=TRUE) library(EpiEstim) } # use estimate_R() funtion, which includes Cori et al (2013)'s method # moving average from https://stackoverflow.com/questions/743812/calculating-moving-average ma <- function(x, n=7, s=1) {filter(x, rep(1 / n, n), sides=s)} x <- covid19(country="JPN", level=2, start="2020-01-01", end=Sys.Date()) y <- subset(x, date=as.Date("2021-01-01")) palette("Okabe-Ito") matplot(as.Date(zz$date, "%d/%m/%y"), cbind(zz$Rt1, zz$Rt2, zz$Rt3, zz$Rt4), xlim=c(min(zz$date), max(zz$date)+20), ylim=c(0, 4), frame=FALSE, col=1:4, lty=1:4, type="l", xlab="時間経過(月)", ylab="Rt", main="東京のCOVID-19について4つの方法で推定したRtの推移 [Data] Guidotti E, Ardia D (2020) COVID-19 Data Hub. Working paper https://doi.org/10.13140/RG.2.2.11649.81763", sub=sprintf("%s まで", max(zz$date))) segments(as.Date("2021-01-01"), 1, max(zz$date)+20, 1, col="red") # segments(as.Date("2020-04-08"), 0, as.Date("2020-04-08"), 4, lty=3, lwd=3, col="darkred") # text(as.Date("2020-04-08"), 4, "緊急事態宣言\n発出", pos=3) # segments(as.Date("2020-05-14"), 0, as.Date("2020-05-14"), 4, lty=3, lwd=3, col="navy") # text(as.Date("2020-05-14"), 4, "39県解除", pos=3) # segments(as.Date("2020-06-19"), 0, as.Date("2020-06-19"), 4, lty=3, lwd=3, col="navy") # text(as.Date("2020-06-19"), 4, "東京解除", pos=3) # segments(as.Date("2020-07-22"), 0, as.Date("2020-07-22"), 4, lty=3, lwd=3, col="navy") # text(as.Date("2020-07-22"), 4, "GoToキャン\nペーン", pos=3) # segments(as.Date("2020-10-01"), 0, as.Date("2020-10-01"), 4, lty=3, lwd=3, col="navy") # text(as.Date("2020-10-01"), 4, "GoTo東京追加\nGoToイート開始", pos=3) segments(as.Date("2021-01-07"), 0, as.Date("2021-01-07"), 4, lty=3, lwd=3, col="darkred") text(as.Date("2021-01-07"), 3, "一都三県緊急\n事態宣言発出", pos=3) segments(as.Date("2021-01-14"), 0, as.Date("2021-01-14"), 4, lty=3, lwd=3, col="darkred") text(as.Date("2021-01-14"), 2, "二府五県緊急\n事態宣言追加", pos=3) segments(as.Date("2021-02-08"), 0, as.Date("2021-02-08"), 4, lty=3, lwd=3, col="navy") text(as.Date("2021-02-08"), 3, "栃木県緊急\n事態宣言解除", pos=3) segments(as.Date("2021-02-28"), 0, as.Date("2021-02-28"), 4, lty=3, lwd=3, col="navy") text(as.Date("2021-02-28"), 3, "二府四県緊急\n事態宣言解除", pos=3) segments(as.Date("2021-03-21"), 0, as.Date("2021-03-21"), 4, lty=3, lwd=3, col="navy") text(as.Date("2021-03-21"), 3, "一都三県緊急\n事態宣言解除", pos=3) segments(as.Date("2021-04-05"), 0, as.Date("2021-04-05"), 4, lty=3, lwd=3, col="darkred") text(as.Date("2021-04-05"), 2, "六市蔓延防止\n重点措置開始", pos=3) segments(as.Date("2021-04-25"), 0, as.Date("2021-04-25"), 4, lty=3, lwd=3, col="darkred") text(as.Date("2021-04-25"), 3, "4都府県緊急\n事態宣言発出", pos=3) segments(as.Date("2021-05-12"), 0, as.Date("2021-05-12"), 4, lty=3, lwd=3, col="darkred") text(as.Date("2021-05-12"), 2, "愛知福岡緊急\n事態宣言発出", pos=3) segments(as.Date("2021-05-16"), 0, as.Date("2021-05-16"), 4, lty=3, lwd=3, col="darkred") text(as.Date("2021-05-16"), 3, "一道二県緊急\n事態宣言発出", pos=3) segments(as.Date("2021-06-20"), 0, as.Date("2021-06-20"), 3, lty=3, lwd=3, col="navy") text(as.Date("2021-06-20"), 2, "沖縄以外緊急\n事態宣言解除", pos=3) segments(as.Date("2021-07-12"), 0, as.Date("2021-07-12"), 4, lty=3, lwd=3, col="darkred") text(as.Date("2021-07-12"), 3, "東京緊急事態\n宣言発出", pos=3) legend("topright", col=1:4, lty=1:4, legend=c("Coriら(2013)の方法", "東洋経済の方法", "感染研の方法", "東洋経済の方法補正無し"))