In kimfinale/pfilterCOVID:
knitr::opts_chunk$set(echo=TRUE, results='hide', message=FALSE, warning=FALSE)
Install the pfilterCOVID package
You will need to install the package, "pfilterCOVID", before using the pfilter function in the document
# devtools::load_all()
# devtools::install()
devtools::install_github("kimfinale/pfilterCOVID", force = TRUE)
Load packages
library(pfilterCOVID)
library(deSolve)
library(tidyverse)
library(cowplot)
library(EpiEstim)
library(foreach)
library(doParallel)
source('R/deconvolve.R')
theme_set(theme_bw())
Set parameter values
# beta is updated by R0
THETA <- list(beta = 0.3, R0 = 1.2, epsilon = 1/2.5, delta = 1/5,
gamma = 1/2.5, bp = 1, ba = 1, fa = 0.3, fd = 1, I0 = 1e2,
time_dep_Rt = FALSE)
# This changes the default y0 in the pfilterCOVID package
THETA$Y0 <- c(S = 1e7 - THETA$I0, E = 0, P = 0, A = 0,
I = THETA$I0, R = 0, CE = 0, CI = 10, CR = 0)
THETA$tend <- 200
usethis::use_data(THETA, overwrite = TRUE)
params <- THETA
Run the ODE version (generate the dataset, D1, and Figure 2A)
params <- THETA
tend <- params$tend
params$I0 <- 100
params$Y0 <- c(S = 1e7 - params$I0, E = 0, P = 0, A = 0,
I = params$I0, R = 0, CE = 0, CI = 0, CR = 0)
tstamp <- format(Sys.time(), "%Y%m%d")
times <- seq(0, tend, by = 1)
params$time_dep_Rt <- TRUE
ode(y = params$Y0, times = times, func = ode_sepair, parms = params) %>%
as.data.frame() -> out
out$daily_infected <- c(0, diff(out$CE))
out$daily_symptom <- c(0, diff(out$CI))
out$daily_confirmed <- c(0, diff(out$CR))
out$daily_Rt <- sapply(0:tend, function(x) get_Rt(x))
# saveRDS(out, paste0("outputs/ode_I0_100_", tstamp, ".rds"))
Figure 2A
tstamp <- "20230619"
out <- readRDS(paste0("outputs/ode_I0_100_", tstamp, ".rds"))
out_long <- out %>% pivot_longer(-time)
out_long$name <- factor(out_long$name,
levels = c("S", "E", "P", "A", "I", "R",
"CE", "CI", "CR", "daily_infected",
"daily_symptom", "daily_confirmed", "daily_Rt"))
fig2a <- out_long %>%
filter(name == "daily_infected" |
name == "daily_symptom" | name == "daily_confirmed") %>%
ggplot(aes(x = time, y = value, color = name)) +
geom_line(linewidth=1.2) +
labs(x='Time (day)', y = 'Number of individuals', color = "") +
theme(legend.position = c(.2, .9),
legend.background = element_rect(fill='transparent')) +
scale_color_hue(labels=c("daily_infected" = "Infection",
"daily_symptom"="Symptom",
"daily_confirmed"="Confirmation"))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
Figure 2B
Plots of Rt and infection/confirmation time series of stoch outputs
NB. Figure 1 is a compartment flow diagram which was not generated by code
# Rt data set; Rt_dataset_I0=100_20220105.rds, Rt_dataset_I0=10_20220105.rds
tstamp <- "20220105"
dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
# Generate Figure 2B.
# plot for the pre-defined Rt
fig2b <- ggplot() +
geom_line(aes(x=dat$time, y = dat$Rt), linewidth=1.4) +
labs(x="Time (day)", y = expression(paste(italic(R)[italic(t)])))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text = element_text(size=13))
figure2 <- plot_grid(fig2a, fig2b, nrow = 2, labels = "AUTO")
figure2
# Generate plots for a comparison of simulated data (deterministic vs. stochastic, both with perfect observation) for daily_infected or confirmed.
oderes <- dat$ode[, c("time", "daily_infected", "daily_confirmed")]
oderes %>% pivot_longer(cols = -time) -> oderes
stochres <- dat$stoch_perfect_obs$daily_confirmed[,1:length(dat$stoch_perfect_obs$daily_confirmed)]
stochres$time <- dat$time
stochres %>% pivot_longer(cols = -time) -> stochres
# plt_case <- ggplot() +
# geom_line(data = stochres, aes(time, value, group=name), color="grey50") +
# geom_line(data = oderes, aes(time, value, color=name), size=1.2) +
# labs(x="time", y="", color="")
Run the Gillespie's direct method (generate the datset, D2, and Figure 3)
This is to create a stochastic data set based on the Gillespie's direct method
set.seed(42)
# tstamp <- format(Sys.time(), "%Y%m%d")
nrun <- 20 # increase the number to obtain large samples
res <- gillespie_run(func = gillespie_sepair,
tend = tend,
nrun = nrun,
y0 = params$Y0,
params = params,
report_dt = 1)
# saveRDS(res, "outputs/res_20230619.rds")
Figure 3
res <- readRDS("outputs/res_20230619.rds")
nrun <- 20
tend <- params$tend
daily_infected <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
daily_symptom <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
daily_confirmed <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
for (i in seq_len(length(res))) {
inf <- diff(res[[i]]$CE)
onset <- diff(res[[i]]$CI)
conf <- diff(res[[i]]$CR)
daily_infected[1:(length(inf)+1), i] <- c(0, inf)
daily_symptom[1:(length(onset)+1), i] <- c(0, onset)
daily_confirmed[1:(length(conf)+1), i] <- c(0, conf)
}
daily_infected <- data.frame(matrix(0, nrow = tend + 1, ncol = nrow(dat$stoch_perfect_obs$daily_infected)))
for (i in 1:nrow(dat$stoch_perfect_obs$daily_infected)) {
daily_infected[,i] <- dat$stoch_perfect_obs$daily_infected[[i]]
}
stoch <- list()
stoch$daily_confirmed <- daily_confirmed
stoch$daily_symptom <- daily_symptom
stoch$daily_infected <- daily_infected
stoch$daily_Rt <- sapply(0:tend, function(x) get_Rt(x))
# "stoch" contains daily time series for infection, symptom, and confirmation (corresponds to the dataset, D2, in the manuscript)
# saveRDS(stoch, paste0("outputs/stoch_I0_100_", tstamp, ".rds"))
df <- daily_infected
df$time <- 0:tend
df %>% pivot_longer(cols = -time) -> df
# Generate Figure 3.
locy = c(2000,1920,1840)
fig3 <- ggplot() +
geom_line(data = df, aes(time, value, group = name),
color = "grey80") +
geom_line(data = out, aes(time, daily_infected), # ODE model output
color = "darkred",
size = 1.5, inherit.aes = FALSE) +
labs(x="Time (day)", y="Daily infected cases", color="") +
geom_segment(aes(x=c(0,0,0), xend=c(5,5,5), y=locy,
yend=locy), size=c(1,1,1),
color = c("grey80", "darkred","black"),
linetype = c("solid","solid","dotted")) +
geom_line(aes(dat$time, rowMeans(dat$stoch_perfect_obs$daily_infected)), # mean of stochastic simulations ## simply to load the data we are using (2,000 simulations) - save time
linetype = "dotted",
size = 1.5, inherit.aes = FALSE) +
geom_text(aes(x=c(7,7,7), y=locy), size=c(4.6,4.6,4.6),
hjust = c(0,0,0), vjust = c(0.5,0.5,0.5),
label = c("Stochastic simulations",
"Determinstic simulation",
"Mean of stochastic simulations"))+
theme(text = element_text(size=16),
axis.text = element_text(size=15),
legend.text=element_text(size=15))
fig3
# ggsave("plots/daily_case_20230619.png", plot=figure2,
# width=3.4*3, height=2.7*3, units="in")
Total dataset
The data set that contains all the necessary information: Rt trajectory, ODE
and stochastic model result
Rt_dataset <- list()
Rt_dataset$README <- "To see if potential bias of Rt estimation based on the confirmation data set, especially in the most recent time, is reduced or eliminated when the number of cases are large"
Rt_dataset$ode <- out
Rt_dataset$stoch_perfect_obs <- stoch
Rt_dataset$Rt <- sapply(0:tend, function(x) get_Rt(x))
Rt_dataset$time <- out$time
Rt_dataset$params <- params
Rt_dataset$y0 <- params$Y0
# saveRDS(Rt_dataset, paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
Rt estimation using the particle filter
Particle filtering - particle_filter
Expect errors since backward sampling is not in place
set.seed(1)
tstamp <- "20220105"
I0 <- 100
dat <- readRDS(paste0("outputs/Rt_dataset_I0=", I0, "_", tstamp, ".rds"))
# determine the data type to model (infection vs. confirmation)
dtype <- "confirmation"
# dtype <- "infection"
if (dtype == "confirmation") {
d <- dat$ode[, c("time", "daily_confirmed")]
} else {
d <- dat$ode[, c("time", "daily_infected")]
}
dat$params$time_dep_Rt <- FALSE # this must be set FALSE to produce daily estimates
params <- dat$params
params[["betavol"]] <- 0.2 #standard deviation of the prior distribution
npart = 1e4 # increase the number of particles to obtain more accurate estimates
pf <- pfilter(
params = params,
y = params$Y0,
data = d,
data_type = dtype,
npart = npart,
tend = nrow(d),
dt = 0.1,
error_pdf = "pois",
systematic_resampling = FALSE,
backward_sampling = TRUE,
stoch = FALSE)
# Rt <- readRDS("outputs/Rt_default.rds")
epsilon <- params[["epsilon"]]
delta <- params[["delta"]]
gamma <- params[["gamma"]]
bp <- params[["bp"]] # relative infectiousness of pre-symptomatic state
ba <- params[["ba"]] # relative infectiousness of asymptomatic state
fa <- params[["fa"]] # fraction of asymptomatic state
durP <- (1 / delta - 1 / epsilon)
durI <- (1 / gamma)
R0_dur <- ((1 - fa) + fa * ba) * durI + bp * durP
med <- apply(pf$beta_filtered, 1, quantile, probs = c(0.5))
# plot(dat$time + 2, dat$daily_Rt, type="l", ylim=c(0, max(upper * R0_dur)))
# Comarison between the pre-determined Rt and PF (without backward sampling) inferred Rt based on ODE daily_confirmed
plot(dat$ode$time + 2, dat$ode$daily_Rt, type="l")
lines(med * R0_dur, lwd=2, col=3)
Fitting the infection or confirmation time series from the ODE model
This uses the extract_trace function that calls pfilter function
and also executes backward sampling
# apply on infection time series
tstamp <- "20220105"
dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
params <- dat$params
params$time_dep_Rt <- FALSE
params[["betavol"]] <- 0.2
# increase the number to improve accuracy
nrep <- 1e3
npart <- 1e4
dt <- 0.2
library(parallel)
library(doParallel)
ncores <- detectCores() - 1
d_inf <- data.frame(date = dat$time, daily_infected = round(dat$ode$daily_infected))
set.seed(42)
cl <- makeCluster(getOption("cl.cores", ncores))
doParallel::registerDoParallel(cl)
PF based on the infection time series
pf_inf <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_inf,
data_type = "infection",
npart = npart,
tend = nrow(d_inf),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_inf, "outputs/pf_inf_20230619.rds")
Figure 4
pf_inf <- readRDS("outputs/pf_inf_20230619.rds")
nrep <- 1e3
npart <- 1e4
dt <- 0.2
dtype <- "infection"
tstamp <- "20220105"
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_inf, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_inf <- cbind(Rt_quantile, d_inf)
col_fill <- "deepskyblue4"
col_dat <- "grey70"
df_inf$Rt <- dat$ode$daily_Rt
fig4a <- ggplot(df_inf[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = col_fill, alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = col_fill, alpha = 0.8) +
geom_line(aes(y = `50%`), color = "darkblue", size = 1.2) +
geom_line(aes(x = date, y = dat$ode$daily_Rt[80:200]), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.25))+scale_y_continuous(breaks=seq(0,2.25,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig4a
# ggsave(paste0("plots/", parset_chr, ".png"), plt)
# apply on confirmation time series
d_conf <- data.frame(date = dat$time,
daily_confirmed = round(dat$ode$daily_confirmed))
PF based on the confirmation time series
set.seed(42)
cl <- makeCluster(getOption("cl.cores", ncores))
doParallel::registerDoParallel(cl)
pf_conf <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_conf,
data_type = "confirmation",
npart = npart,
tend = nrow(d_conf),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_conf, "outputs/pf_conf_20230619.rds")
Figure 4b
pf_conf <- readRDS("outputs/pf_conf_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
Rt_est <- as.data.frame(sapply(pf_conf, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_conf <- cbind(Rt_quantile, d_conf)
df_conf$Rt <- dat$ode$daily_Rt
fig4b <- ggplot(df_conf[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "orange", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "brown", size = 1.2) +
geom_line(aes(x = date, y = dat$ode$daily_Rt[80:200]), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.25))+scale_y_continuous(breaks=seq(0,2.25,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig4b
# put all together for figure 3
figure4 <- plot_grid(fig4a, fig4b, nrow = 2, labels = "AUTO")
figure4
# ggsave("plots/Rt_recovered_20230619.png", plot=figure3,
# width=3.4*2, height=2.7*2, units="in")
Fitting confirmation time series from the stochastic model
tstamp <- "20220105"
dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
params <- dat$params
params$time_dep_Rt <- FALSE
params[["betavol"]] <- 0.2
nrep <- 1e3
npart <- 1e4
dt <- 0.2
d_stoc_19 <- data.frame(date = dat$time, daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,19]))
library(parallel)
library(doParallel)
set.seed(42)
ncores <- detectCores()
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
# apply on the stochastic seed #19 (fig 4b)
pf_stoc_19 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_19,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_19),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_19, "outputs/pf_stoc_19_20230619.rds")
Figure 5b
pf_stoc_19 <- readRDS("outputs/pf_stoc_19_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_19, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_19 <- cbind(Rt_quantile, d_stoc_19)
col_fill <- "red"
col_dat <- "grey70"
df_stoc_19$Rt <- dat$ode$daily_Rt
fig5b <- ggplot(df_stoc_19[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = col_fill, alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = col_fill, alpha = 0.8) +
geom_line(aes(y = `50%`), color = "darkred", size = 1.2) +
geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+
scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5b
PF based on the confirmation time series derived from the stochastic model
# apply on the stochastic seed #04 (fig 4c)
d_stoc_04 <- data.frame(date = dat$time,
daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,4]))
set.seed(42)
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
pf_stoc_04 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_04,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_04),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_04, "outputs/pf_stoc_04_20230619.rds")
pf_stoc_04 <- readRDS("outputs/pf_stoc_04_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_04, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_04 <- cbind(Rt_quantile, d_stoc_04)
col_fill <- "green"
col_dat <- "grey70"
df_stoc_04$Rt <- dat$ode$daily_Rt
fig5c <- ggplot(df_stoc_04[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "#76ff7a", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "green", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "darkgreen", size = 1.2) +
geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.5))+
scale_y_continuous(breaks=seq(0,2.5,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5c
PF based on the confirmation time series dervied from the stochastic model
Incidence level is lower than the deterministic model
# apply on the stochastic seed #02 (fig 4d)
d_stoc_02 <- data.frame(date = dat$time,
daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,2]))
set.seed(42)
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
pf_stoc_02 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_02,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_02),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_02, "outputs/pf_stoc_02_20230619.rds")
Figure 5d
pf_stoc_02 <- readRDS("outputs/pf_stoc_02_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_02, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_02 <- cbind(Rt_quantile, d_stoc_02)
col_fill <- "blue"
col_dat <- "grey70"
df_stoc_02$Rt <- dat$ode$daily_Rt
fig5d <- ggplot(df_stoc_02[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "orange", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "brown", size = 1.2) +
geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.5))+
scale_y_continuous(breaks=seq(0, 2.5, by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5d
# Generate figure 5
df_stoc <- data.frame(date=dat$ode$time, ode=dat$ode$daily_confirmed, c2=dat$stoch_perfect_obs$daily_confirmed[,2], c4=dat$stoch_perfect_obs$daily_confirmed[,4], c19=dat$stoch_perfect_obs$daily_confirmed[,19])
fig5a <- ggplot(df_stoc) + geom_line(aes(x=date, y=ode), linewidth=1, color="black")+geom_line(aes(x=date,y=c2),linewidth=1,color="orange")+geom_line(aes(x=date,y=c4),color="green")+geom_line(aes(x=date,y=c19),color="red") +
labs(y = "Number of confirmation", x = "Time (day)") +
scale_color_hue(labels=c("ode" = "ODE", "c2"="Low", "c4"="Medium", "c19"="High")) +
geom_segment(aes(x=0, xend=5, y=1000, yend=1000),
size=1.1, color = "red", linetype = "solid") +
geom_segment(aes(x=0, xend=5, y=950, yend=950),
size=1.1, color = "green", linetype = "solid") +
geom_segment(aes(x=0, xend=5, y=900, yend=900),
size=1.1, color = "orange", linetype = "solid") +
geom_text(aes(x=7, y=1000), size=4.6, hjust=0, vjust=0.5, label="High intensity") +
geom_text(aes(x=7, y=950), size=4.6, hjust=0, vjust=0.5, label="Medium intensity") +
geom_text(aes(x=7, y=900), size=4.6, hjust=0, vjust=0.5, label="Low intensity")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5a
left_col <- plot_grid(fig5a, labels = c("A"))
right_col <- plot_grid(fig5b, fig5c, fig5d, nrow = 3, labels = c("B", "C", "D"))
figure5 <- plot_grid(left_col, right_col, ncol = 2)
figure5
# ggsave("plots/Rt_stoch_recovered_20230619.png", plot=figure4,
# width=3.4*2, height=2.7*2, units="in")
Rt estimation using EpiEstim
Computing mean and std of generation time
# infection (disease) age for generation time
norm_factor <- integrate(function(x) {
1 - pgamma(x, shape = 2, rate = 1/2.5)}, 0, Inf)
f <- function(x) {
(1 - pgamma(x, shape = 2, rate = 1/2.5))/norm_factor$value}
# latent time
g <- function(x) { dexp(x, rate = 1/2.5) }
z <- function(z) {
integrate(function(x, z) g(z-x)*f(x), 0, Inf, z)$value}
fz <- Vectorize(z)
## upper is set to arbitrarily large number
mean_int <- integrate(function(y) y*fz(y), lower = 0,
upper = 100, subdivisions = 1e4)
var_int <- integrate(function(y) y^2*fz(y), lower = 0,
upper = 100, subdivisions = 1e4)
var <- var_int$value - mean_int$value^2 # variance = E[x^2] - E[X]^2
mean_int$value
sqrt(var)
# from Mathematica
mean <- 6.25 # mean_si
std <- 4.14578 # std_si
Apply EpiEstim on daily confirmation series
Expect discrepancy (i.e., delay) between the pre-defined Rt and the estimated Rt)
Figure 5
# Moving average
MA <- function(input_ts, lookback) {
sum <- 0
look_back1 <- lookback + 1
ma_input <- rep(0, length(input_ts))
for (i in 1:length(input_ts)) {
sum = sum + input_ts[i]
if (i < look_back1) {
ma_input[i] <- sum / i
} else {
sum = sum - input_ts[i-7]
ma_input[i] <- sum / 7
}
}
return(ma_input)
}
# apply EpiEstim on confirmation time series (ODE version)
ma_input <- MA(dat$ode$daily_confirmed, 7)
T <- length(ma_input)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(ma_input, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
df_epi_ode_conf <- data.frame(sapply(res_daily$R,c))
df_epi_ode_conf$Rt <- dat$ode$daily_Rt[3:201]
df_epi_ode_conf$time <- t_start-1
col_fill <- "deepskyblue4"
col_dat <- "grey70"
figure6 <-
ggplot(df_epi_ode_conf[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkblue", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure6
# ggsave("plots/EpiEstim_Rt_delay_20230619.png", plot=figure5,
# width=3.4*2, height=2.7, units="in")
Introduce a delay distribution to address the delay observed above (for deconvolutions)
delay_dist <- {
t <- seq(0, 19)
p_unnorm <- pgamma(t + 1, shape = 3, rate = 1/2.5) -
pgamma(t, shape = 3, rate = 1/2.5)
data.frame(delay = t, pmf = p_unnorm / sum(p_unnorm))
}
Rt based on renewal equation (without misspecification effect of SI)
Figure 7
Estimate Rt by applying deconvolution on confimration series, followed by EpiEstim
# apply on the stochastic seed #19 (fig 6a)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,19]
# obtain moving average before deconvolution
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
# output to the above opearation is defined as "result$x_unobs"
unobs <- result$x_unobs
# simple plots for the output of deconvolution (i.e., inferred infection), and the infection from the ODE model
# determine an appropriate lag value by comparing the peaks
# plot(unobs)
# lines(dat$ode$daily_infected)
lag <- which.max(dat$ode$daily_infected) - which.max(unobs)
#lag
# adjust the inferred infection (i.e., unobs)
# exclude first 4 elements of unobs (when, e.g., lag=4) as out of estimation period of our interest
# this causes 4 less estimates (i.e., upto t=196 instead of t=200 for lag=4)
unobs <- unobs[5:length(unobs)]
# now apply EpiEstim
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
# used the mean and std of generation time above for mean_si and std_si
df_epi_19 <- data.frame(sapply(res_daily$R,c))
df_epi_19 <- rbind(df_epi_19, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_19$time <- dat$time[3:201]
df_epi_19$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$sto_per_high <- df_epi_19$Median.R.
plt_epi_19 <-
ggplot(df_epi_19[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#ff6961", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "red", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkred", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.5, 2))+
scale_y_continuous(breaks=seq(0,2,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #4 (fig 6b)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
# graph
df_epi_4 <- data.frame(sapply(res_daily$R,c))
df_epi_4 <- rbind(df_epi_4, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4$time <- dat$time[3:201]
df_epi_4$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$sto_per_high <- df_epi_4$Median.R.
plt_epi_4 <-
ggplot(df_epi_4[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#76ff7a", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "green", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkgreen", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.5, 2))+
scale_y_continuous(breaks=seq(0,2,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #2 (fig 6c)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,2]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
df_epi_2 <- data.frame(sapply(res_daily$R,c))
df_epi_2 <- rbind(df_epi_2, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_2$time <- dat$time[3:201]
df_epi_2$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$sto_per_high <- df_epi_2$Median.R.
plt_epi_2 <-
ggplot(df_epi_2[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "brown", size = 1) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "darkblue", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.5, 2))+
scale_y_continuous(breaks=seq(0,2,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure7 <- plot_grid(plt_epi_19, plt_epi_4, plt_epi_2, nrow = 3, labels = c("A", "B", "C"))
figure7
# ggsave("plots/EpiEstim_Rt_deconvolution_20230619.png", plot=figure6,
# width=3.4*2, height=2.7*2, units="in")
Rt based on renewal equation (with misspecification effect of SI)
Figure 8
Estimate Rt by applying deconvolution on confimration series, followed by EpiEstim
# apply on the stochastic seed #19 (fig 7a)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,19]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
# misspecification effect: used mean_si=5, istead of 6.25
df_epi_19m <- data.frame(sapply(res_daily$R,c))
df_epi_19m <- rbind(df_epi_19m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_19m$time <- dat$time[3:201]
df_epi_19m$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$mis_sto_per_low <- df_epi_2m$Median.R.
plt_epi_19m <-
ggplot(df_epi_19m[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#ff6961", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "red", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkred", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #4 (fig 7b)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
# output to the above opearation is defined as "result$x_unobs"
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
# epiestim_4m
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <-
estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
df_epi_4m <- data.frame(sapply(res_daily$R,c))
df_epi_4m <- rbind(df_epi_4m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4m$time <- dat$time[3:201]
df_epi_4m$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$mis_sto_per_low <- df_epi_2m$Median.R.
plt_epi_4m <-
ggplot(df_epi_4m[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#76ff7a", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "green", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkgreen", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #2 (fig 7c)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,2]
ma <- seq(0, length(input_ts)-1)
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
df_epi_2m <- data.frame(sapply(res_daily$R,c))
df_epi_2m <- rbind(df_epi_2m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_2m$time <- dat$time[3:201]
df_epi_2m$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$mis_sto_per_low <- df_epi_2m$Median.R.
plt_epi_2m <-
ggplot(df_epi_2m[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "brown", size = 1) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "darkblue", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure8 <- plot_grid(plt_epi_19m, plt_epi_4m, plt_epi_2m, nrow = 3, labels = c("a", "b", "c"))
figure8
# ggsave("plots/EpiEstim_stoch_Rt_deconvolution_20230619.png", plot=figure7,
# width=3.4*2, height=2.7*2, units="in")
Performance comarison between partifle filtering and renewal equation-based estimations
Figure 9
Plot for the curve of pre-defined Rt by periods of different epidemic characteristics
R_flat <- c(dat$Rt[1:101], rep(NA,32), dat$Rt[134:195], rep(NA, 6))
R_curve <- c(rep(NA, 100), dat$Rt[101:127], rep(NA, 74))
R_curve_last <- c(rep(NA, 126), dat$Rt[127:134], rep(NA, 67))
R_last <- c(rep(NA, 194), dat$Rt[195:201])
figure9 <- ggplot()+
geom_line(aes(x=dat$time, y=R_flat), color="black", linewidth=1.2)+
geom_line(aes(x=dat$time, y=R_curve), color="red", linewidth=1.2)+
geom_line(aes(x=dat$time, y=R_curve_last), linetype="dotted", color="red", linewidth=1.2)+
geom_line(aes(x=dat$time, y=R_last), linetype="dotted", linewidth=1.2)+
labs(y=expression(italic(R)[italic(t)]), x="Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure9
# ggsave("plots/Rt_sections_20230619.png", plot=figure8,
# width=3.4*2, height=2.7*2, units="in")
Extra simulaions to terminate estimation on the last day of the dynamic period (i.e., day = 132)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4][1:133]
# PF
d_stoc_04_mid <- data.frame(date = dat$time[1:133],
daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,4][1:133]))
set.seed(42)
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
pf_stoc_04_mid <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_04_mid,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_04_mid),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_04_mid, "outputs/pf_stoc_04_mid_20230619.rds")
Figure 10
pf_stoc_04_mid <- readRDS("outputs/pf_stoc_04_mid_20230619.rds")
parset_chr <- paste0("confirmation", "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_04_mid, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_04_mid <- cbind(Rt_quantile, d_stoc_04_mid)
col_fill <- "deepskyblue4"
col_dat <- "grey70"
df_stoc_04_mid$Rt <- dat$ode$daily_Rt[1:133]
fig10a <- ggplot(df_stoc_04_mid[81:133,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill="orange", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill="orange", alpha = 0.8) +
geom_line(aes(y = `50%`), color="brown", size = 1.2) +
geom_line(aes(x = date, y = df_stoc_04_mid$Rt[80:132]), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# EpiEstim
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time[1:133]),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
# adjustment of unobs
unobs <- unobs[5:133]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
# res_daily: without misspecification of SI
# res_daily_m: with misspecification of SI
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
res_daily_m <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
df_epi_4_mid <- data.frame(sapply(res_daily$R, c))
df_epi_4_mid_m <- data.frame(sapply(res_daily_m$R, c))
df_epi_4_mid <- df_epi_4_mid[,c(5,7,8,9,11)]
df_epi_4_mid_m <- df_epi_4_mid_m[,c(5,7,8,9,11)]
df_epi_4_mid <- rbind(df_epi_4_mid, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4_mid_m <- rbind(df_epi_4_mid_m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4_mid$time <- dat$time[3:133]
df_epi_4_mid$Rt <- dat$ode$daily_Rt[3:133]
df_epi_4_mid_m$time <- dat$time[3:133]
df_epi_4_mid_m$Rt <- dat$ode$daily_Rt[3:133]
# fig10b
fig10b <-
ggplot(df_epi_4_mid[79:131,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = col_fill, alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = col_fill, alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkblue", size = 1.2) +
geom_line(aes(y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+
scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig10c
fig10c <-
ggplot(df_epi_4_mid_m[79:131,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.),
fill="orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax=Quantile.0.75.R.),
fill="orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "brown", size = 1.2) +
geom_line(aes(y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+
scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure10 <- plot_grid(fig10a, fig10b, fig10c, nrow=3, labels=c("A", "B", "C"))
figure10
# ggsave("plots/Rt_PF_EpiEstim_20230619.png", plot=figure9,
# width=3.4*2, height=2.7*2, units="in")
Calculate RMSE between the pre-defined Rt and estimated Rt's (for Table 2)
dg <- df_stoc_04[82:201,]
dg1 <- df_epi_4[84:199,]
dg2 <- df_epi_4m[84:199,]
# #
flat_Rt <- c(dg$Rt[1:20], dg$Rt[53:120])
flat_PF_Rt <- c(dg$`50%`[1:20], dg$`50%`[53:120])
flat_epi_Rt <- c(dg1$Median.R.[1:20], dg1$Median.R.[53:116], NA, NA, NA, NA)
flat_epim_Rt <- c(dg2$Median.R.[1:20], dg2$Median.R.[53:116], NA, NA, NA, NA)
# #
general_period <- list("PF" = sqrt(mean((df_stoc_04[82:201,]$Rt - df_stoc_04[82:201,]$`50%`)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((df_stoc_04[82:197,]$Rt - df_epi_4[84:199,]$Median.R.)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((df_stoc_04[82:197,]$Rt - df_epi_4m[84:199,]$Median.R.)^2, na.rm=TRUE)))
flat_period <- list("PF" = sqrt(mean((flat_Rt - flat_PF_Rt)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((flat_Rt - flat_epi_Rt)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((flat_Rt - flat_epim_Rt)^2, na.rm=TRUE)))
last_flat <- list("PF" = sqrt(mean((flat_Rt[82:88] - flat_PF_Rt[82:88])^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((flat_Rt[82:88] - flat_epi_Rt[82:88])^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((flat_Rt[82:88] - flat_epim_Rt[82:88])^2, na.rm=TRUE)))
fluctuation <- list("PF" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_stoc_04_mid[99:133,]$`50%`)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_epi_4_mid[99:133,]$Median.R.)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_epi_4_mid_m[99:133,]$Median.R.)^2, na.rm=TRUE)))
last_fluctuaion <- list("PF" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_stoc_04_mid[127:133,]$`50%`)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_epi_4_mid[127:133,]$Median.R.)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_epi_4_mid_m[127:133,]$Median.R.)^2, na.rm=TRUE)))
# RMSE table (in the order of PF, EpiEstim,and EpiEstim with misspecification)
table2 <- list("general period" = general_period, "flat period" = flat_period, "last 1 week flat" = last_flat, "fluctuation period" = fluctuation, "last 1 week fluctuation" = last_fluctuaion)
Rt for COVID-19 in Korea during 1 Sep 2020 - 31 Oct 2021
set.seed(42)
cl <- makeCluster(getOption("cl.cores", detectCores()))
registerDoParallel(cl)
nrep <- 1000
npart <- 10000
params[["betavol"]] <- 0.20
y0 <- params$Y0
y0["I"] <- 100
y0["E"] <- 100
data <- read_csv("data/daily_confirmed_202009_202110.csv")
# data <- data[order(data$date),]
# write_csv(data, "data/daily_confirmed_202009_202110.csv")
# data$date <- 0:(nrow(data)-1)
pf_conf <- foreach(i=1:nrep, .packages=c("pfilterCOVID"), .inorder=F) %dopar% {
extract_trace(params = params,
y = y0,
data = data,
data_type = "confirmation",
npart = npart,
tend = nrow(data),
dt = 0.1,
error_pdf = "negbin",
negbin_size = 100,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_conf, "outputs/pf_conf_covid_korea_20230729.rds")
Generate Figure 11 (supplmentary figure)
Red and green vertical lines represent strengthening and weakening of the social distancing measures.
Figure 11
dtype <- "confirmation"
tstamp <- format(Sys.Date(), "%Y%m%d")
data <- read_csv("data/daily_confirmed_202009_202110.csv")
pf_korea <- readRDS("outputs/pf_conf_covid_korea_20230729.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_korea, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
cumul_est <- as.data.frame(sapply(pf_korea, function(x) x[, "CR"]))
cumul_quantile <- as.data.frame(t(apply(cumul_est, 1, function(x) quantile(x, pr))))
Rtdf <- cbind(Rt_quantile, data)
Rtdf$date <- data$date
cumuldf <- cbind(cumul_quantile, data)
cumuldf$date <- data$date
npi_data <- data.frame(date=as.Date(c("2020-10-11", "2021-03-12", "2021-06-11",
"2020-11-17", "2020-11-22", "2020-12-05", "2020-12-08", "2020-12-21")),
degree = c("Low","Low","Low","High","High","High","High","High"))
npi_high <- npi_data %>% filter(degree=="High")
npi_low <- npi_data %>% filter(degree=="Low")
# Rt
fig11a <- ggplot(Rtdf, aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "steelblue", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "steelblue", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "steelblue", size = 1.2) +
geom_hline(yintercept = 1, color = "black", size = 1, linetype="dotted")+
geom_vline(xintercept=npi_high$date, color="firebrick", size=1)+
geom_vline(xintercept=npi_low$date, color="darkgreen", size=1)+
labs(y = expression(italic(R)[italic(t)]), x="") +
scale_x_date(limits=as.Date(c("2020-10-01","2021-06-30")), date_breaks = "2 week")+
scale_y_continuous(limits=c(0,3),breaks=seq(0,3,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13),
axis.text.x=element_text(angle=90,hjust=1))
# fig11a
fig11b <- ggplot(cumuldf) +
geom_col(aes(x=date, y=daily_confirmed), fill="grey50", alpha=0.3) +
geom_ribbon(aes(x=date, ymin=`2.5%`, ymax=`97.5%`),
fill="steelblue", alpha=0.3) +
geom_ribbon(aes(x = date, ymin = `25%`, ymax=`75%`),
fill="steelblue", alpha=0.7) +
geom_line(aes(x=date, y=`50%`), color="steelblue", size=1)+
# scale_y_continuous(limits=c(0,2000),breaks=seq(0,1000,by=200))+
labs(title = "", x = "", y = "Daily confirmed case") +
scale_x_date(limits=as.Date(c("2020-10-01","2021-06-30")), date_breaks = "2 week")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13),
axis.text.x=element_text(angle=90,hjust=1))
# fig11b
figure11 <- plot_grid(fig11a, fig11b, nrow=2, labels="AUTO")
figure11
# ggsave("plots/Rt_confimed_case_korea_20230730.png", plot=figure11,
# width=3.4*2, height=2.7*3, units="in")
kimfinale/pfilterCOVID documentation built on July 30, 2023, 12:15 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo=TRUE, results='hide', message=FALSE, warning=FALSE)
Install the pfilterCOVID package
You will need to install the package, "pfilterCOVID", before using the pfilter function in the document
# devtools::load_all() # devtools::install() devtools::install_github("kimfinale/pfilterCOVID", force = TRUE)
Load packages
library(pfilterCOVID) library(deSolve) library(tidyverse) library(cowplot) library(EpiEstim) library(foreach) library(doParallel) source('R/deconvolve.R') theme_set(theme_bw())
Set parameter values
# beta is updated by R0 THETA <- list(beta = 0.3, R0 = 1.2, epsilon = 1/2.5, delta = 1/5, gamma = 1/2.5, bp = 1, ba = 1, fa = 0.3, fd = 1, I0 = 1e2, time_dep_Rt = FALSE) # This changes the default y0 in the pfilterCOVID package THETA$Y0 <- c(S = 1e7 - THETA$I0, E = 0, P = 0, A = 0, I = THETA$I0, R = 0, CE = 0, CI = 10, CR = 0) THETA$tend <- 200 usethis::use_data(THETA, overwrite = TRUE) params <- THETA
Run the ODE version (generate the dataset, D1, and Figure 2A)
params <- THETA tend <- params$tend params$I0 <- 100 params$Y0 <- c(S = 1e7 - params$I0, E = 0, P = 0, A = 0, I = params$I0, R = 0, CE = 0, CI = 0, CR = 0) tstamp <- format(Sys.time(), "%Y%m%d") times <- seq(0, tend, by = 1) params$time_dep_Rt <- TRUE ode(y = params$Y0, times = times, func = ode_sepair, parms = params) %>% as.data.frame() -> out out$daily_infected <- c(0, diff(out$CE)) out$daily_symptom <- c(0, diff(out$CI)) out$daily_confirmed <- c(0, diff(out$CR)) out$daily_Rt <- sapply(0:tend, function(x) get_Rt(x)) # saveRDS(out, paste0("outputs/ode_I0_100_", tstamp, ".rds"))
Figure 2A
tstamp <- "20230619" out <- readRDS(paste0("outputs/ode_I0_100_", tstamp, ".rds")) out_long <- out %>% pivot_longer(-time) out_long$name <- factor(out_long$name, levels = c("S", "E", "P", "A", "I", "R", "CE", "CI", "CR", "daily_infected", "daily_symptom", "daily_confirmed", "daily_Rt")) fig2a <- out_long %>% filter(name == "daily_infected" | name == "daily_symptom" | name == "daily_confirmed") %>% ggplot(aes(x = time, y = value, color = name)) + geom_line(linewidth=1.2) + labs(x='Time (day)', y = 'Number of individuals', color = "") + theme(legend.position = c(.2, .9), legend.background = element_rect(fill='transparent')) + scale_color_hue(labels=c("daily_infected" = "Infection", "daily_symptom"="Symptom", "daily_confirmed"="Confirmation"))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13))
Figure 2B
Plots of Rt and infection/confirmation time series of stoch outputs NB. Figure 1 is a compartment flow diagram which was not generated by code
# Rt data set; Rt_dataset_I0=100_20220105.rds, Rt_dataset_I0=10_20220105.rds tstamp <- "20220105" dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds")) # Generate Figure 2B. # plot for the pre-defined Rt fig2b <- ggplot() + geom_line(aes(x=dat$time, y = dat$Rt), linewidth=1.4) + labs(x="Time (day)", y = expression(paste(italic(R)[italic(t)])))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text = element_text(size=13)) figure2 <- plot_grid(fig2a, fig2b, nrow = 2, labels = "AUTO") figure2 # Generate plots for a comparison of simulated data (deterministic vs. stochastic, both with perfect observation) for daily_infected or confirmed. oderes <- dat$ode[, c("time", "daily_infected", "daily_confirmed")] oderes %>% pivot_longer(cols = -time) -> oderes stochres <- dat$stoch_perfect_obs$daily_confirmed[,1:length(dat$stoch_perfect_obs$daily_confirmed)] stochres$time <- dat$time stochres %>% pivot_longer(cols = -time) -> stochres # plt_case <- ggplot() + # geom_line(data = stochres, aes(time, value, group=name), color="grey50") + # geom_line(data = oderes, aes(time, value, color=name), size=1.2) + # labs(x="time", y="", color="")
Run the Gillespie's direct method (generate the datset, D2, and Figure 3)
This is to create a stochastic data set based on the Gillespie's direct method
set.seed(42) # tstamp <- format(Sys.time(), "%Y%m%d") nrun <- 20 # increase the number to obtain large samples res <- gillespie_run(func = gillespie_sepair, tend = tend, nrun = nrun, y0 = params$Y0, params = params, report_dt = 1) # saveRDS(res, "outputs/res_20230619.rds")
Figure 3
res <- readRDS("outputs/res_20230619.rds") nrun <- 20 tend <- params$tend daily_infected <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun)) daily_symptom <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun)) daily_confirmed <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun)) for (i in seq_len(length(res))) { inf <- diff(res[[i]]$CE) onset <- diff(res[[i]]$CI) conf <- diff(res[[i]]$CR) daily_infected[1:(length(inf)+1), i] <- c(0, inf) daily_symptom[1:(length(onset)+1), i] <- c(0, onset) daily_confirmed[1:(length(conf)+1), i] <- c(0, conf) } daily_infected <- data.frame(matrix(0, nrow = tend + 1, ncol = nrow(dat$stoch_perfect_obs$daily_infected))) for (i in 1:nrow(dat$stoch_perfect_obs$daily_infected)) { daily_infected[,i] <- dat$stoch_perfect_obs$daily_infected[[i]] } stoch <- list() stoch$daily_confirmed <- daily_confirmed stoch$daily_symptom <- daily_symptom stoch$daily_infected <- daily_infected stoch$daily_Rt <- sapply(0:tend, function(x) get_Rt(x)) # "stoch" contains daily time series for infection, symptom, and confirmation (corresponds to the dataset, D2, in the manuscript) # saveRDS(stoch, paste0("outputs/stoch_I0_100_", tstamp, ".rds")) df <- daily_infected df$time <- 0:tend df %>% pivot_longer(cols = -time) -> df # Generate Figure 3. locy = c(2000,1920,1840) fig3 <- ggplot() + geom_line(data = df, aes(time, value, group = name), color = "grey80") + geom_line(data = out, aes(time, daily_infected), # ODE model output color = "darkred", size = 1.5, inherit.aes = FALSE) + labs(x="Time (day)", y="Daily infected cases", color="") + geom_segment(aes(x=c(0,0,0), xend=c(5,5,5), y=locy, yend=locy), size=c(1,1,1), color = c("grey80", "darkred","black"), linetype = c("solid","solid","dotted")) + geom_line(aes(dat$time, rowMeans(dat$stoch_perfect_obs$daily_infected)), # mean of stochastic simulations ## simply to load the data we are using (2,000 simulations) - save time linetype = "dotted", size = 1.5, inherit.aes = FALSE) + geom_text(aes(x=c(7,7,7), y=locy), size=c(4.6,4.6,4.6), hjust = c(0,0,0), vjust = c(0.5,0.5,0.5), label = c("Stochastic simulations", "Determinstic simulation", "Mean of stochastic simulations"))+ theme(text = element_text(size=16), axis.text = element_text(size=15), legend.text=element_text(size=15)) fig3 # ggsave("plots/daily_case_20230619.png", plot=figure2, # width=3.4*3, height=2.7*3, units="in")
Total dataset
The data set that contains all the necessary information: Rt trajectory, ODE and stochastic model result
Rt_dataset <- list() Rt_dataset$README <- "To see if potential bias of Rt estimation based on the confirmation data set, especially in the most recent time, is reduced or eliminated when the number of cases are large" Rt_dataset$ode <- out Rt_dataset$stoch_perfect_obs <- stoch Rt_dataset$Rt <- sapply(0:tend, function(x) get_Rt(x)) Rt_dataset$time <- out$time Rt_dataset$params <- params Rt_dataset$y0 <- params$Y0 # saveRDS(Rt_dataset, paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
Rt estimation using the particle filter
Particle filtering - particle_filter Expect errors since backward sampling is not in place
set.seed(1) tstamp <- "20220105" I0 <- 100 dat <- readRDS(paste0("outputs/Rt_dataset_I0=", I0, "_", tstamp, ".rds")) # determine the data type to model (infection vs. confirmation) dtype <- "confirmation" # dtype <- "infection" if (dtype == "confirmation") { d <- dat$ode[, c("time", "daily_confirmed")] } else { d <- dat$ode[, c("time", "daily_infected")] } dat$params$time_dep_Rt <- FALSE # this must be set FALSE to produce daily estimates params <- dat$params params[["betavol"]] <- 0.2 #standard deviation of the prior distribution npart = 1e4 # increase the number of particles to obtain more accurate estimates pf <- pfilter( params = params, y = params$Y0, data = d, data_type = dtype, npart = npart, tend = nrow(d), dt = 0.1, error_pdf = "pois", systematic_resampling = FALSE, backward_sampling = TRUE, stoch = FALSE) # Rt <- readRDS("outputs/Rt_default.rds") epsilon <- params[["epsilon"]] delta <- params[["delta"]] gamma <- params[["gamma"]] bp <- params[["bp"]] # relative infectiousness of pre-symptomatic state ba <- params[["ba"]] # relative infectiousness of asymptomatic state fa <- params[["fa"]] # fraction of asymptomatic state durP <- (1 / delta - 1 / epsilon) durI <- (1 / gamma) R0_dur <- ((1 - fa) + fa * ba) * durI + bp * durP med <- apply(pf$beta_filtered, 1, quantile, probs = c(0.5)) # plot(dat$time + 2, dat$daily_Rt, type="l", ylim=c(0, max(upper * R0_dur))) # Comarison between the pre-determined Rt and PF (without backward sampling) inferred Rt based on ODE daily_confirmed plot(dat$ode$time + 2, dat$ode$daily_Rt, type="l") lines(med * R0_dur, lwd=2, col=3)
Fitting the infection or confirmation time series from the ODE model
This uses the extract_trace function that calls pfilter function and also executes backward sampling
# apply on infection time series tstamp <- "20220105" dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds")) params <- dat$params params$time_dep_Rt <- FALSE params[["betavol"]] <- 0.2 # increase the number to improve accuracy nrep <- 1e3 npart <- 1e4 dt <- 0.2 library(parallel) library(doParallel) ncores <- detectCores() - 1 d_inf <- data.frame(date = dat$time, daily_infected = round(dat$ode$daily_infected)) set.seed(42) cl <- makeCluster(getOption("cl.cores", ncores)) doParallel::registerDoParallel(cl)
PF based on the infection time series
pf_inf <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% { extract_trace(params = params, y = params$Y0, data = d_inf, data_type = "infection", npart = npart, tend = nrow(d_inf), dt = dt, error_pdf = "pois", negbin_size = 15, backward_sampling = TRUE, stoch = FALSE) } parallel::stopCluster(cl) # saveRDS(pf_inf, "outputs/pf_inf_20230619.rds")
Figure 4
pf_inf <- readRDS("outputs/pf_inf_20230619.rds") nrep <- 1e3 npart <- 1e4 dt <- 0.2 dtype <- "infection" tstamp <- "20220105" parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart, "_nrep=", nrep, "_dt=", dt, "_", tstamp) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_est <- as.data.frame(sapply(pf_inf, function(x) x[, "Rt"])) Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr)))) df_inf <- cbind(Rt_quantile, d_inf) col_fill <- "deepskyblue4" col_dat <- "grey70" df_inf$Rt <- dat$ode$daily_Rt fig4a <- ggplot(df_inf[81:201,], aes(x = date)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = col_fill, alpha = 0.5) + geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = col_fill, alpha = 0.8) + geom_line(aes(y = `50%`), color = "darkblue", size = 1.2) + geom_line(aes(x = date, y = dat$ode$daily_Rt[80:200]), size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.25))+scale_y_continuous(breaks=seq(0,2.25,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # fig4a # ggsave(paste0("plots/", parset_chr, ".png"), plt)
# apply on confirmation time series d_conf <- data.frame(date = dat$time, daily_confirmed = round(dat$ode$daily_confirmed))
PF based on the confirmation time series
set.seed(42) cl <- makeCluster(getOption("cl.cores", ncores)) doParallel::registerDoParallel(cl) pf_conf <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% { extract_trace(params = params, y = params$Y0, data = d_conf, data_type = "confirmation", npart = npart, tend = nrow(d_conf), dt = dt, error_pdf = "pois", negbin_size = 15, backward_sampling = TRUE, stoch = FALSE) } parallel::stopCluster(cl) # saveRDS(pf_conf, "outputs/pf_conf_20230619.rds")
Figure 4b
pf_conf <- readRDS("outputs/pf_conf_20230619.rds") parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart, "_nrep=", nrep, "_dt=", dt, "_", tstamp) Rt_est <- as.data.frame(sapply(pf_conf, function(x) x[, "Rt"])) Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr)))) df_conf <- cbind(Rt_quantile, d_conf) df_conf$Rt <- dat$ode$daily_Rt fig4b <- ggplot(df_conf[81:201,], aes(x = date)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "orange", alpha = 0.5) + geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "orange", alpha = 0.8) + geom_line(aes(y = `50%`), color = "brown", size = 1.2) + geom_line(aes(x = date, y = dat$ode$daily_Rt[80:200]), size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.25))+scale_y_continuous(breaks=seq(0,2.25,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # fig4b # put all together for figure 3 figure4 <- plot_grid(fig4a, fig4b, nrow = 2, labels = "AUTO") figure4 # ggsave("plots/Rt_recovered_20230619.png", plot=figure3, # width=3.4*2, height=2.7*2, units="in")
Fitting confirmation time series from the stochastic model
tstamp <- "20220105" dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds")) params <- dat$params params$time_dep_Rt <- FALSE params[["betavol"]] <- 0.2 nrep <- 1e3 npart <- 1e4 dt <- 0.2 d_stoc_19 <- data.frame(date = dat$time, daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,19]))
library(parallel) library(doParallel) set.seed(42) ncores <- detectCores() cl <- makeCluster(getOption("cl.cores", 2)) doParallel::registerDoParallel(cl) # apply on the stochastic seed #19 (fig 4b) pf_stoc_19 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% { extract_trace(params = params, y = params$Y0, data = d_stoc_19, data_type = "confirmation", npart = npart, tend = nrow(d_stoc_19), dt = dt, error_pdf = "pois", negbin_size = 15, backward_sampling = TRUE, stoch = FALSE) } parallel::stopCluster(cl) # saveRDS(pf_stoc_19, "outputs/pf_stoc_19_20230619.rds")
Figure 5b
pf_stoc_19 <- readRDS("outputs/pf_stoc_19_20230619.rds") parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart, "_nrep=", nrep, "_dt=", dt, "_", tstamp) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_est <- as.data.frame(sapply(pf_stoc_19, function(x) x[, "Rt"])) Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr)))) df_stoc_19 <- cbind(Rt_quantile, d_stoc_19) col_fill <- "red" col_dat <- "grey70" df_stoc_19$Rt <- dat$ode$daily_Rt fig5b <- ggplot(df_stoc_19[81:201,], aes(x = date)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = col_fill, alpha = 0.5) + geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = col_fill, alpha = 0.8) + geom_line(aes(y = `50%`), color = "darkred", size = 1.2) + geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.75))+ scale_y_continuous(breaks=seq(0,2.75,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # fig5b
PF based on the confirmation time series derived from the stochastic model
# apply on the stochastic seed #04 (fig 4c) d_stoc_04 <- data.frame(date = dat$time, daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,4]))
set.seed(42) cl <- makeCluster(getOption("cl.cores", 2)) doParallel::registerDoParallel(cl) pf_stoc_04 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% { extract_trace(params = params, y = params$Y0, data = d_stoc_04, data_type = "confirmation", npart = npart, tend = nrow(d_stoc_04), dt = dt, error_pdf = "pois", negbin_size = 15, backward_sampling = TRUE, stoch = FALSE) } parallel::stopCluster(cl) # saveRDS(pf_stoc_04, "outputs/pf_stoc_04_20230619.rds")
pf_stoc_04 <- readRDS("outputs/pf_stoc_04_20230619.rds") parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart, "_nrep=", nrep, "_dt=", dt, "_", tstamp) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_est <- as.data.frame(sapply(pf_stoc_04, function(x) x[, "Rt"])) Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr)))) df_stoc_04 <- cbind(Rt_quantile, d_stoc_04) col_fill <- "green" col_dat <- "grey70" df_stoc_04$Rt <- dat$ode$daily_Rt fig5c <- ggplot(df_stoc_04[81:201,], aes(x = date)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "#76ff7a", alpha = 0.5) + geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "green", alpha = 0.8) + geom_line(aes(y = `50%`), color = "darkgreen", size = 1.2) + geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.5))+ scale_y_continuous(breaks=seq(0,2.5,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # fig5c
PF based on the confirmation time series dervied from the stochastic model
Incidence level is lower than the deterministic model
# apply on the stochastic seed #02 (fig 4d) d_stoc_02 <- data.frame(date = dat$time, daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,2]))
set.seed(42) cl <- makeCluster(getOption("cl.cores", 2)) doParallel::registerDoParallel(cl) pf_stoc_02 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% { extract_trace(params = params, y = params$Y0, data = d_stoc_02, data_type = "confirmation", npart = npart, tend = nrow(d_stoc_02), dt = dt, error_pdf = "pois", negbin_size = 15, backward_sampling = TRUE, stoch = FALSE) } parallel::stopCluster(cl) # saveRDS(pf_stoc_02, "outputs/pf_stoc_02_20230619.rds")
Figure 5d
pf_stoc_02 <- readRDS("outputs/pf_stoc_02_20230619.rds") parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart, "_nrep=", nrep, "_dt=", dt, "_", tstamp) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_est <- as.data.frame(sapply(pf_stoc_02, function(x) x[, "Rt"])) Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr)))) df_stoc_02 <- cbind(Rt_quantile, d_stoc_02) col_fill <- "blue" col_dat <- "grey70" df_stoc_02$Rt <- dat$ode$daily_Rt fig5d <- ggplot(df_stoc_02[81:201,], aes(x = date)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "orange", alpha = 0.5) + geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "orange", alpha = 0.8) + geom_line(aes(y = `50%`), color = "brown", size = 1.2) + geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.5))+ scale_y_continuous(breaks=seq(0, 2.5, by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # fig5d # Generate figure 5 df_stoc <- data.frame(date=dat$ode$time, ode=dat$ode$daily_confirmed, c2=dat$stoch_perfect_obs$daily_confirmed[,2], c4=dat$stoch_perfect_obs$daily_confirmed[,4], c19=dat$stoch_perfect_obs$daily_confirmed[,19]) fig5a <- ggplot(df_stoc) + geom_line(aes(x=date, y=ode), linewidth=1, color="black")+geom_line(aes(x=date,y=c2),linewidth=1,color="orange")+geom_line(aes(x=date,y=c4),color="green")+geom_line(aes(x=date,y=c19),color="red") + labs(y = "Number of confirmation", x = "Time (day)") + scale_color_hue(labels=c("ode" = "ODE", "c2"="Low", "c4"="Medium", "c19"="High")) + geom_segment(aes(x=0, xend=5, y=1000, yend=1000), size=1.1, color = "red", linetype = "solid") + geom_segment(aes(x=0, xend=5, y=950, yend=950), size=1.1, color = "green", linetype = "solid") + geom_segment(aes(x=0, xend=5, y=900, yend=900), size=1.1, color = "orange", linetype = "solid") + geom_text(aes(x=7, y=1000), size=4.6, hjust=0, vjust=0.5, label="High intensity") + geom_text(aes(x=7, y=950), size=4.6, hjust=0, vjust=0.5, label="Medium intensity") + geom_text(aes(x=7, y=900), size=4.6, hjust=0, vjust=0.5, label="Low intensity")+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # fig5a left_col <- plot_grid(fig5a, labels = c("A")) right_col <- plot_grid(fig5b, fig5c, fig5d, nrow = 3, labels = c("B", "C", "D")) figure5 <- plot_grid(left_col, right_col, ncol = 2) figure5 # ggsave("plots/Rt_stoch_recovered_20230619.png", plot=figure4, # width=3.4*2, height=2.7*2, units="in")
Rt estimation using EpiEstim
Computing mean and std of generation time
# infection (disease) age for generation time norm_factor <- integrate(function(x) { 1 - pgamma(x, shape = 2, rate = 1/2.5)}, 0, Inf) f <- function(x) { (1 - pgamma(x, shape = 2, rate = 1/2.5))/norm_factor$value} # latent time g <- function(x) { dexp(x, rate = 1/2.5) } z <- function(z) { integrate(function(x, z) g(z-x)*f(x), 0, Inf, z)$value} fz <- Vectorize(z) ## upper is set to arbitrarily large number mean_int <- integrate(function(y) y*fz(y), lower = 0, upper = 100, subdivisions = 1e4) var_int <- integrate(function(y) y^2*fz(y), lower = 0, upper = 100, subdivisions = 1e4) var <- var_int$value - mean_int$value^2 # variance = E[x^2] - E[X]^2 mean_int$value sqrt(var) # from Mathematica mean <- 6.25 # mean_si std <- 4.14578 # std_si
Apply EpiEstim on daily confirmation series
Expect discrepancy (i.e., delay) between the pre-defined Rt and the estimated Rt)
Figure 5
# Moving average MA <- function(input_ts, lookback) { sum <- 0 look_back1 <- lookback + 1 ma_input <- rep(0, length(input_ts)) for (i in 1:length(input_ts)) { sum = sum + input_ts[i] if (i < look_back1) { ma_input[i] <- sum / i } else { sum = sum - input_ts[i-7] ma_input[i] <- sum / 7 } } return(ma_input) } # apply EpiEstim on confirmation time series (ODE version) ma_input <- MA(dat$ode$daily_confirmed, 7) T <- length(ma_input) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 res_daily <- estimate_R(ma_input, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578)))) df_epi_ode_conf <- data.frame(sapply(res_daily$R,c)) df_epi_ode_conf$Rt <- dat$ode$daily_Rt[3:201] df_epi_ode_conf$time <- t_start-1 col_fill <- "deepskyblue4" col_dat <- "grey70" figure6 <- ggplot(df_epi_ode_conf[80:199,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "darkblue", size = 1.2) + geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) figure6 # ggsave("plots/EpiEstim_Rt_delay_20230619.png", plot=figure5, # width=3.4*2, height=2.7, units="in")
Introduce a delay distribution to address the delay observed above (for deconvolutions)
delay_dist <- { t <- seq(0, 19) p_unnorm <- pgamma(t + 1, shape = 3, rate = 1/2.5) - pgamma(t, shape = 3, rate = 1/2.5) data.frame(delay = t, pmf = p_unnorm / sum(p_unnorm)) }
Rt based on renewal equation (without misspecification effect of SI)
Figure 7
Estimate Rt by applying deconvolution on confimration series, followed by EpiEstim
# apply on the stochastic seed #19 (fig 6a) input_ts <- dat$stoch_perfect_obs$daily_confirmed[,19] # obtain moving average before deconvolution ma <- seq(0, length(input_ts)-1) lookback <- 7 sum <- 0 for (i in 1:length(input_ts)){ sum = sum + input_ts[i] if (i < lookback+1) { ma[i] = sum / i } else { sum = sum - input_ts[i-7] ma[i] = sum / 7 } } result <- deconvolve_rltype_goldsteinetal( t_obs_min = dat$time[1], y_obs = round(ma), delay_min = delay_dist$delay[1], pmf_delay = delay_dist$pmf, t_unobs_min = dat$time[1], n_unobs = length(dat$time), n_iterations_max = 100, n_iterations = NULL ) # output to the above opearation is defined as "result$x_unobs" unobs <- result$x_unobs # simple plots for the output of deconvolution (i.e., inferred infection), and the infection from the ODE model # determine an appropriate lag value by comparing the peaks # plot(unobs) # lines(dat$ode$daily_infected) lag <- which.max(dat$ode$daily_infected) - which.max(unobs) #lag # adjust the inferred infection (i.e., unobs) # exclude first 4 elements of unobs (when, e.g., lag=4) as out of estimation period of our interest # this causes 4 less estimates (i.e., upto t=196 instead of t=200 for lag=4) unobs <- unobs[5:length(unobs)] # now apply EpiEstim T <- length(unobs) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 res_daily <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578)))) # used the mean and std of generation time above for mean_si and std_si df_epi_19 <- data.frame(sapply(res_daily$R,c)) df_epi_19 <- rbind(df_epi_19, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_19$time <- dat$time[3:201] df_epi_19$Rt <- dat$ode$daily_Rt[3:201] # df_deconv$sto_per_high <- df_epi_19$Median.R. plt_epi_19 <- ggplot(df_epi_19[80:199,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#ff6961", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "red", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "darkred", size = 1.2) + geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.5, 2))+ scale_y_continuous(breaks=seq(0,2,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # apply on the stochastic seed #4 (fig 6b) input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4] ma <- seq(0, length(input_ts)-1) lookback <- 7 sum <- 0 for (i in 1:length(input_ts)){ sum = sum + input_ts[i] if (i < lookback+1) { ma[i] = sum / i } else { sum = sum - input_ts[i-7] ma[i] = sum / 7 } } result <- deconvolve_rltype_goldsteinetal( t_obs_min = dat$time[1], y_obs = round(ma), delay_min = delay_dist$delay[1], pmf_delay = delay_dist$pmf, t_unobs_min = dat$time[1], n_unobs = length(dat$time), n_iterations_max = 100, n_iterations = NULL ) unobs <- result$x_unobs unobs <- unobs[5:length(unobs)] T <- length(unobs) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 res_daily <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578)))) # graph df_epi_4 <- data.frame(sapply(res_daily$R,c)) df_epi_4 <- rbind(df_epi_4, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_4$time <- dat$time[3:201] df_epi_4$Rt <- dat$ode$daily_Rt[3:201] # df_deconv$sto_per_high <- df_epi_4$Median.R. plt_epi_4 <- ggplot(df_epi_4[80:199,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#76ff7a", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "green", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "darkgreen", size = 1.2) + geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.5, 2))+ scale_y_continuous(breaks=seq(0,2,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # apply on the stochastic seed #2 (fig 6c) input_ts <- dat$stoch_perfect_obs$daily_confirmed[,2] ma <- seq(0, length(input_ts)-1) lookback <- 7 sum <- 0 for (i in 1:length(input_ts)){ sum = sum + input_ts[i] if (i < lookback+1) { ma[i] = sum / i } else { sum = sum - input_ts[i-7] ma[i] = sum / 7 } } result <- deconvolve_rltype_goldsteinetal( t_obs_min = dat$time[1], y_obs = round(ma), delay_min = delay_dist$delay[1], pmf_delay = delay_dist$pmf, t_unobs_min = dat$time[1], n_unobs = length(dat$time), n_iterations_max = 100, n_iterations = NULL ) unobs <- result$x_unobs unobs <- unobs[5:length(unobs)] T <- length(unobs) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 res_daily <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578)))) df_epi_2 <- data.frame(sapply(res_daily$R,c)) df_epi_2 <- rbind(df_epi_2, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_2$time <- dat$time[3:201] df_epi_2$Rt <- dat$ode$daily_Rt[3:201] # df_deconv$sto_per_high <- df_epi_2$Median.R. plt_epi_2 <- ggplot(df_epi_2[80:199,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "brown", size = 1) + geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "darkblue", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.5, 2))+ scale_y_continuous(breaks=seq(0,2,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) figure7 <- plot_grid(plt_epi_19, plt_epi_4, plt_epi_2, nrow = 3, labels = c("A", "B", "C")) figure7 # ggsave("plots/EpiEstim_Rt_deconvolution_20230619.png", plot=figure6, # width=3.4*2, height=2.7*2, units="in")
Rt based on renewal equation (with misspecification effect of SI)
Figure 8
Estimate Rt by applying deconvolution on confimration series, followed by EpiEstim
# apply on the stochastic seed #19 (fig 7a) input_ts <- dat$stoch_perfect_obs$daily_confirmed[,19] ma <- seq(0, length(input_ts)-1) lookback <- 7 sum <- 0 for (i in 1:length(input_ts)){ sum = sum + input_ts[i] if (i < lookback+1) { ma[i] = sum / i } else { sum = sum - input_ts[i-7] ma[i] = sum / 7 } } result <- deconvolve_rltype_goldsteinetal( t_obs_min = dat$time[1], y_obs = round(ma), delay_min = delay_dist$delay[1], pmf_delay = delay_dist$pmf, t_unobs_min = dat$time[1], n_unobs = length(dat$time), n_iterations_max = 100, n_iterations = NULL ) unobs <- result$x_unobs unobs <- unobs[5:length(unobs)] T <- length(unobs) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 res_daily <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578)))) # misspecification effect: used mean_si=5, istead of 6.25 df_epi_19m <- data.frame(sapply(res_daily$R,c)) df_epi_19m <- rbind(df_epi_19m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_19m$time <- dat$time[3:201] df_epi_19m$Rt <- dat$ode$daily_Rt[3:201] # df_deconv$mis_sto_per_low <- df_epi_2m$Median.R. plt_epi_19m <- ggplot(df_epi_19m[80:199,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#ff6961", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "red", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "darkred", size = 1.2) + geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # apply on the stochastic seed #4 (fig 7b) input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4] ma <- seq(0, length(input_ts)-1) lookback <- 7 sum <- 0 for (i in 1:length(input_ts)){ sum = sum + input_ts[i] if (i < lookback+1) { ma[i] = sum / i } else { sum = sum - input_ts[i-7] ma[i] = sum / 7 } } result <- deconvolve_rltype_goldsteinetal( t_obs_min = dat$time[1], y_obs = round(ma), delay_min = delay_dist$delay[1], pmf_delay = delay_dist$pmf, t_unobs_min = dat$time[1], n_unobs = length(dat$time), n_iterations_max = 100, n_iterations = NULL ) # output to the above opearation is defined as "result$x_unobs" unobs <- result$x_unobs unobs <- unobs[5:length(unobs)] # epiestim_4m T <- length(unobs) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 res_daily <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578)))) df_epi_4m <- data.frame(sapply(res_daily$R,c)) df_epi_4m <- rbind(df_epi_4m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_4m$time <- dat$time[3:201] df_epi_4m$Rt <- dat$ode$daily_Rt[3:201] # df_deconv$mis_sto_per_low <- df_epi_2m$Median.R. plt_epi_4m <- ggplot(df_epi_4m[80:199,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#76ff7a", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "green", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "darkgreen", size = 1.2) + geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # apply on the stochastic seed #2 (fig 7c) input_ts <- dat$stoch_perfect_obs$daily_confirmed[,2] ma <- seq(0, length(input_ts)-1) sum <- 0 for (i in 1:length(input_ts)){ sum = sum + input_ts[i] if (i < lookback+1) { ma[i] = sum / i } else { sum = sum - input_ts[i-7] ma[i] = sum / 7 } } result <- deconvolve_rltype_goldsteinetal( t_obs_min = dat$time[1], y_obs = round(ma), delay_min = delay_dist$delay[1], pmf_delay = delay_dist$pmf, t_unobs_min = dat$time[1], n_unobs = length(dat$time), n_iterations_max = 100, n_iterations = NULL ) unobs <- result$x_unobs unobs <- unobs[5:length(unobs)] T <- length(unobs) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 res_daily <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578)))) df_epi_2m <- data.frame(sapply(res_daily$R,c)) df_epi_2m <- rbind(df_epi_2m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_2m$time <- dat$time[3:201] df_epi_2m$Rt <- dat$ode$daily_Rt[3:201] # df_deconv$mis_sto_per_low <- df_epi_2m$Median.R. plt_epi_2m <- ggplot(df_epi_2m[80:199,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "brown", size = 1) + geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "darkblue", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) figure8 <- plot_grid(plt_epi_19m, plt_epi_4m, plt_epi_2m, nrow = 3, labels = c("a", "b", "c")) figure8 # ggsave("plots/EpiEstim_stoch_Rt_deconvolution_20230619.png", plot=figure7, # width=3.4*2, height=2.7*2, units="in")
Performance comarison between partifle filtering and renewal equation-based estimations
Figure 9
Plot for the curve of pre-defined Rt by periods of different epidemic characteristics
R_flat <- c(dat$Rt[1:101], rep(NA,32), dat$Rt[134:195], rep(NA, 6)) R_curve <- c(rep(NA, 100), dat$Rt[101:127], rep(NA, 74)) R_curve_last <- c(rep(NA, 126), dat$Rt[127:134], rep(NA, 67)) R_last <- c(rep(NA, 194), dat$Rt[195:201]) figure9 <- ggplot()+ geom_line(aes(x=dat$time, y=R_flat), color="black", linewidth=1.2)+ geom_line(aes(x=dat$time, y=R_curve), color="red", linewidth=1.2)+ geom_line(aes(x=dat$time, y=R_curve_last), linetype="dotted", color="red", linewidth=1.2)+ geom_line(aes(x=dat$time, y=R_last), linetype="dotted", linewidth=1.2)+ labs(y=expression(italic(R)[italic(t)]), x="Time (day)")+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) figure9 # ggsave("plots/Rt_sections_20230619.png", plot=figure8, # width=3.4*2, height=2.7*2, units="in")
Extra simulaions to terminate estimation on the last day of the dynamic period (i.e., day = 132)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4][1:133] # PF d_stoc_04_mid <- data.frame(date = dat$time[1:133], daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,4][1:133]))
set.seed(42) cl <- makeCluster(getOption("cl.cores", 2)) doParallel::registerDoParallel(cl) pf_stoc_04_mid <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% { extract_trace(params = params, y = params$Y0, data = d_stoc_04_mid, data_type = "confirmation", npart = npart, tend = nrow(d_stoc_04_mid), dt = dt, error_pdf = "pois", negbin_size = 15, backward_sampling = TRUE, stoch = FALSE) } parallel::stopCluster(cl) # saveRDS(pf_stoc_04_mid, "outputs/pf_stoc_04_mid_20230619.rds")
Figure 10
pf_stoc_04_mid <- readRDS("outputs/pf_stoc_04_mid_20230619.rds") parset_chr <- paste0("confirmation", "_I0=", params$I0, "_npart=", npart, "_nrep=", nrep, "_dt=", dt, "_", tstamp) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_est <- as.data.frame(sapply(pf_stoc_04_mid, function(x) x[, "Rt"])) Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr)))) df_stoc_04_mid <- cbind(Rt_quantile, d_stoc_04_mid) col_fill <- "deepskyblue4" col_dat <- "grey70" df_stoc_04_mid$Rt <- dat$ode$daily_Rt[1:133] fig10a <- ggplot(df_stoc_04_mid[81:133,], aes(x = date)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill="orange", alpha = 0.5) + geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill="orange", alpha = 0.8) + geom_line(aes(y = `50%`), color="brown", size = 1.2) + geom_line(aes(x = date, y = df_stoc_04_mid$Rt[80:132]), size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.75))+scale_y_continuous(breaks=seq(0,2.75,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # EpiEstim ma <- seq(0, length(input_ts)-1) lookback <- 7 sum <- 0 for (i in 1:length(input_ts)){ sum = sum + input_ts[i] if (i < lookback+1) { ma[i] = sum / i } else { sum = sum - input_ts[i-7] ma[i] = sum / 7 } } result <- deconvolve_rltype_goldsteinetal( t_obs_min = dat$time[1], y_obs = round(ma), delay_min = delay_dist$delay[1], pmf_delay = delay_dist$pmf, t_unobs_min = dat$time[1], n_unobs = length(dat$time[1:133]), n_iterations_max = 100, n_iterations = NULL ) unobs <- result$x_unobs # adjustment of unobs unobs <- unobs[5:133] T <- length(unobs) time_window <- 1 t_start <- seq(3, T-(time_window-1)) t_end <- t_start + time_window-1 # res_daily: without misspecification of SI # res_daily_m: with misspecification of SI res_daily <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578)))) res_daily_m <- estimate_R(unobs, method = "parametric_si", config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578)))) df_epi_4_mid <- data.frame(sapply(res_daily$R, c)) df_epi_4_mid_m <- data.frame(sapply(res_daily_m$R, c)) df_epi_4_mid <- df_epi_4_mid[,c(5,7,8,9,11)] df_epi_4_mid_m <- df_epi_4_mid_m[,c(5,7,8,9,11)] df_epi_4_mid <- rbind(df_epi_4_mid, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_4_mid_m <- rbind(df_epi_4_mid_m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA)) df_epi_4_mid$time <- dat$time[3:133] df_epi_4_mid$Rt <- dat$ode$daily_Rt[3:133] df_epi_4_mid_m$time <- dat$time[3:133] df_epi_4_mid_m$Rt <- dat$ode$daily_Rt[3:133] # fig10b fig10b <- ggplot(df_epi_4_mid[79:131,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = col_fill, alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = col_fill, alpha = 0.8) + geom_line(aes(y = Median.R.), color = "darkblue", size = 1.2) + geom_line(aes(y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.75))+ scale_y_continuous(breaks=seq(0,2.75,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) # fig10c fig10c <- ggplot(df_epi_4_mid_m[79:131,], aes(x = time)) + geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill="orange", alpha = 0.5) + geom_ribbon(aes(ymin = Quantile.0.25.R., ymax=Quantile.0.75.R.), fill="orange", alpha = 0.8) + geom_line(aes(y = Median.R.), color = "brown", size = 1.2) + geom_line(aes(y = Rt), colour = "black", size = 1, linetype = "dashed") + geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+ labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") + coord_cartesian(ylim = c(0.25, 2.75))+ scale_y_continuous(breaks=seq(0,2.75,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13)) figure10 <- plot_grid(fig10a, fig10b, fig10c, nrow=3, labels=c("A", "B", "C")) figure10 # ggsave("plots/Rt_PF_EpiEstim_20230619.png", plot=figure9, # width=3.4*2, height=2.7*2, units="in")
Calculate RMSE between the pre-defined Rt and estimated Rt's (for Table 2)
dg <- df_stoc_04[82:201,] dg1 <- df_epi_4[84:199,] dg2 <- df_epi_4m[84:199,] # # flat_Rt <- c(dg$Rt[1:20], dg$Rt[53:120]) flat_PF_Rt <- c(dg$`50%`[1:20], dg$`50%`[53:120]) flat_epi_Rt <- c(dg1$Median.R.[1:20], dg1$Median.R.[53:116], NA, NA, NA, NA) flat_epim_Rt <- c(dg2$Median.R.[1:20], dg2$Median.R.[53:116], NA, NA, NA, NA) # # general_period <- list("PF" = sqrt(mean((df_stoc_04[82:201,]$Rt - df_stoc_04[82:201,]$`50%`)^2, na.rm=TRUE)), "EpiEstim" = sqrt(mean((df_stoc_04[82:197,]$Rt - df_epi_4[84:199,]$Median.R.)^2, na.rm=TRUE)), "EpiEstim_m" = sqrt(mean((df_stoc_04[82:197,]$Rt - df_epi_4m[84:199,]$Median.R.)^2, na.rm=TRUE))) flat_period <- list("PF" = sqrt(mean((flat_Rt - flat_PF_Rt)^2, na.rm=TRUE)), "EpiEstim" = sqrt(mean((flat_Rt - flat_epi_Rt)^2, na.rm=TRUE)), "EpiEstim_m" = sqrt(mean((flat_Rt - flat_epim_Rt)^2, na.rm=TRUE))) last_flat <- list("PF" = sqrt(mean((flat_Rt[82:88] - flat_PF_Rt[82:88])^2, na.rm=TRUE)), "EpiEstim" = sqrt(mean((flat_Rt[82:88] - flat_epi_Rt[82:88])^2, na.rm=TRUE)), "EpiEstim_m" = sqrt(mean((flat_Rt[82:88] - flat_epim_Rt[82:88])^2, na.rm=TRUE))) fluctuation <- list("PF" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_stoc_04_mid[99:133,]$`50%`)^2, na.rm=TRUE)), "EpiEstim" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_epi_4_mid[99:133,]$Median.R.)^2, na.rm=TRUE)), "EpiEstim_m" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_epi_4_mid_m[99:133,]$Median.R.)^2, na.rm=TRUE))) last_fluctuaion <- list("PF" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_stoc_04_mid[127:133,]$`50%`)^2, na.rm=TRUE)), "EpiEstim" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_epi_4_mid[127:133,]$Median.R.)^2, na.rm=TRUE)), "EpiEstim_m" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_epi_4_mid_m[127:133,]$Median.R.)^2, na.rm=TRUE))) # RMSE table (in the order of PF, EpiEstim,and EpiEstim with misspecification) table2 <- list("general period" = general_period, "flat period" = flat_period, "last 1 week flat" = last_flat, "fluctuation period" = fluctuation, "last 1 week fluctuation" = last_fluctuaion)
Rt for COVID-19 in Korea during 1 Sep 2020 - 31 Oct 2021
set.seed(42) cl <- makeCluster(getOption("cl.cores", detectCores())) registerDoParallel(cl) nrep <- 1000 npart <- 10000 params[["betavol"]] <- 0.20 y0 <- params$Y0 y0["I"] <- 100 y0["E"] <- 100 data <- read_csv("data/daily_confirmed_202009_202110.csv") # data <- data[order(data$date),] # write_csv(data, "data/daily_confirmed_202009_202110.csv") # data$date <- 0:(nrow(data)-1) pf_conf <- foreach(i=1:nrep, .packages=c("pfilterCOVID"), .inorder=F) %dopar% { extract_trace(params = params, y = y0, data = data, data_type = "confirmation", npart = npart, tend = nrow(data), dt = 0.1, error_pdf = "negbin", negbin_size = 100, backward_sampling = TRUE, stoch = FALSE) } parallel::stopCluster(cl) # saveRDS(pf_conf, "outputs/pf_conf_covid_korea_20230729.rds")
Generate Figure 11 (supplmentary figure) Red and green vertical lines represent strengthening and weakening of the social distancing measures.
Figure 11
dtype <- "confirmation" tstamp <- format(Sys.Date(), "%Y%m%d") data <- read_csv("data/daily_confirmed_202009_202110.csv") pf_korea <- readRDS("outputs/pf_conf_covid_korea_20230729.rds") parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart, "_nrep=", nrep, "_dt=", dt, "_", tstamp) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_est <- as.data.frame(sapply(pf_korea, function(x) x[, "Rt"])) Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr)))) cumul_est <- as.data.frame(sapply(pf_korea, function(x) x[, "CR"])) cumul_quantile <- as.data.frame(t(apply(cumul_est, 1, function(x) quantile(x, pr)))) Rtdf <- cbind(Rt_quantile, data) Rtdf$date <- data$date cumuldf <- cbind(cumul_quantile, data) cumuldf$date <- data$date npi_data <- data.frame(date=as.Date(c("2020-10-11", "2021-03-12", "2021-06-11", "2020-11-17", "2020-11-22", "2020-12-05", "2020-12-08", "2020-12-21")), degree = c("Low","Low","Low","High","High","High","High","High")) npi_high <- npi_data %>% filter(degree=="High") npi_low <- npi_data %>% filter(degree=="Low") # Rt fig11a <- ggplot(Rtdf, aes(x = date)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "steelblue", alpha = 0.5) + geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "steelblue", alpha = 0.8) + geom_line(aes(y = `50%`), color = "steelblue", size = 1.2) + geom_hline(yintercept = 1, color = "black", size = 1, linetype="dotted")+ geom_vline(xintercept=npi_high$date, color="firebrick", size=1)+ geom_vline(xintercept=npi_low$date, color="darkgreen", size=1)+ labs(y = expression(italic(R)[italic(t)]), x="") + scale_x_date(limits=as.Date(c("2020-10-01","2021-06-30")), date_breaks = "2 week")+ scale_y_continuous(limits=c(0,3),breaks=seq(0,3,by=0.5))+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13), axis.text.x=element_text(angle=90,hjust=1)) # fig11a fig11b <- ggplot(cumuldf) + geom_col(aes(x=date, y=daily_confirmed), fill="grey50", alpha=0.3) + geom_ribbon(aes(x=date, ymin=`2.5%`, ymax=`97.5%`), fill="steelblue", alpha=0.3) + geom_ribbon(aes(x = date, ymin = `25%`, ymax=`75%`), fill="steelblue", alpha=0.7) + geom_line(aes(x=date, y=`50%`), color="steelblue", size=1)+ # scale_y_continuous(limits=c(0,2000),breaks=seq(0,1000,by=200))+ labs(title = "", x = "", y = "Daily confirmed case") + scale_x_date(limits=as.Date(c("2020-10-01","2021-06-30")), date_breaks = "2 week")+ theme(text = element_text(size=16), axis.text = element_text(size=13), legend.text=element_text(size=13), axis.text.x=element_text(angle=90,hjust=1)) # fig11b figure11 <- plot_grid(fig11a, fig11b, nrow=2, labels="AUTO") figure11 # ggsave("plots/Rt_confimed_case_korea_20230730.png", plot=figure11, # width=3.4*2, height=2.7*3, units="in")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.